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PART I: INTRODUCTION 


1 . 1 General 

Among various choices involved in applying the finite element 
method for the resolution of problems of fluid or solid mechanics, 
the choice of the element type associated with a given mesh and that 
of the integration rule used on these elements often decide the pro- 
perties of the approximation of the operator to be discretized. A 
bad approximation generally results in a kernel larger than the ker- 
nel of the continuous operator, and these spurious elements of this 
kernel may appear in the discrete solution. This solution then can 
exhibit undesirable oscillations and instabilities. 

In particular, for the general class of problem: 

Au - B*p = f 
Bu = 0 

oscillations and instabilities may arise from a bad approximation 
of the governing operator A or of the constraint operator B. 

As far as the constraint is conerned, numerous theoretical and 
numerical studies of pressure instabilities have been done over the 
last decade [7, 10, 33, 36, 45]. In these studies, it has been pro- 
ven that a key stability condition, the L.B.B. condition, must be 
satisfied in order to have existence and uniqueness of a solution 
[31, 9, 2], and that the constant appearing in this condition must 
be independent of the mesh size in order to have stability of the 
element considered [39, 40, 41]. Several stable and unstable elements 
have been studied with respect to the L.B.B. condition [19, 20, 42,43, 
44]. 
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Our purpose in this report is to concentrate on a bad finite 
element approximation of the governing operator obtained when under- 
integration is used in numerical code for several model problems: 
the Poisson problem, the linear elasticity problem, and for problems 
in the nonlinear theory of elasticity. For each of these problems, 
the reasons for the occurrence of instabilities will be given, way 
to control or eliminate them will be presented, and theorems of exis- 
tence, uniqueness and convergence for the given methods will be estab- 
lished. Finally, numerical results are given which illustrate the 
theory . 

1 . 2 Major Results 

Historical background as well as precise definitions and notations 
will be given in the introductions of Parts II and III. However, we 
list here the major results we have obtained: 

1) In underintegrated finite element methods, a rank-deficiency 
counting for each element cannot work to predict the exact number of 
spurious modes of the gobal stiffness matrix. 

2) When the spurious modes are known exactly, it is posssible 

to eliminate them a-posteriori and obtain a unique, stable and accurate 
solution. 

3) Even when the spurious zero-energy modes cannot be predicted 
from rank-deficiency observations, oscillations can be predicted from 
the eigenvalue analysis of the stiffness matrix. 

4) Concentrated forces are the most effective source of oscilla- 
tions; a way to handle such loads without exciting oscillations is 
theoretically obtained and its numerical implementation turns out to 


be efficient. 
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3) In linear elasticity, the spurious modes and their construc- 
tion is described. 

6) Stabilization methods proposed by several authors may not 

work. 

7) For a general mesh, an element-by-element spurious mode con- 
trol is presented, as well as its numerical results: this control is 

efficient and leads to an accurate solution. 

8) This control is cheap, easy to implement, and can be used 
independently of the material. 

9) In nonlinear incompressible elasticity, several numerical 
results are presented and confirm that the results previously listed 
can be extended to nonlinear problems. 

1.3 Report Organization 

Part II is devoted to the analysis and control of hourglass modes 
in underintegrated finite element methods. The framework of the theo- 
retical study is the simple two-dimensional Poisson problem solved 
using bilinear elements; the comparison involves the 4- and 1-integration 
point rules. We prove that for some boundary conditions, the stiffness 
matrix is rank-deficient and spurious modes appear, but we also prove 
that the solution obtained from the underintegrated problem can be 
processed a-posteriori in order to obtain a stable and convergent solu- 
tion in which the spurious modes have been eliminated. Moreover, the 
method of proof allows us to demonstrate precisely why spurious modes 
are excited, even though they cannot be predicted by a rank— deficiency 
of the stiffness matrix. Numerical experiments are discussed which 
not only illustrate the theory, but also generalize the results to 
various elements (8- and 9-node elements) and operators (linear elasticity 


operator) . 



In Part III, we numerically investigate the behavior of under- 
integrated 8- and 9-node elements associated with linear discontinuous 
pressures for the analysis of problems in finite elasticity. We ob- 
serve that whereas the 9-node element is stable, the 8-node element 
exhibits pressure oscillations. We study the performance of the con- 
trol presented in Part II to the control of oscillations appearing in 
the finite element solution of rubber elasticity when underintegration 


is practiced. 



PART II: ANALYSIS OF INSTABILITIES IN 


UNDERINTEGRATED FINITE ELEMENT METHODS 
2.1 Introduction 

For many years, a special type of numerical instability has been 
observed in finite difference approximations of flow fields, which has 
been referred to as "hourglassing", "keystoning", or "chickenwiring". 
These graphic terms refer to geometrical patterns which appear in com- 
puted flow fields (e.g. velocities) and which emerge as spurious oscil- 
lations superimposed on an otherwise smooth field, the spurious oscilla- 
tions often taking a zig-zag form which resembles an hourglass or a 
chickenwire mesh. These spurious modes can be amplified upon refining 
the mesh, and to control such numerical instabilities, various schemes 
for incorporating "hourglass viscosity" or "hourglass damping" have 
been proposed by some authors. 

It is now known that hourglass modes can arise from an incomplete 
(or poor) approximation of the kernel of the operators in the momentum 
equations in flow or solid mechanics problems (or, more generally, of 
the principal part of the operator in the governing differential equa- 
tion of a given boundary-value problem). For example, in addition to 
the rigid body motions residing in the kernel of the standard operators 
appearing in the equilibrium (momentum) equations of solid and fluid 
mechanics, one finds hourglass modes in various crude discrete models 
of these operators. 

In recent years, the occurrence of hourglass instabilities in 
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under integrated finite element approximations has been observed. In 
the implementation of most finite element methods, integrals defining 
stiffness matrices are evaluated using numerical quadrature schemes. 

To improve computational efficiency, the practice of uTiderintegvcition 
is often employed, by which is meant the use of a quadrature rule of 
an order lower than that required to integrate polynomial integrands 
exactly. This can produce rank-deficient stiffness matrices or, equi- 
valently, an expanded kernel of the equilibrium operation which contains 
spurious hourglass modes, and the result is again a numerically un- 
stable scheme. 

In order to overcome this difficulty, artificial stiffness or vis- 
cosity methods, or other stabilization methods have been proposed by 
several authors (e.g. [3-6, 18, 30 ] ) . These methods involve computing 
an under integrated matrix, and then adding a stabilization matrix which 
effectively eliminates the hourglass modes. They turn out to be fairly 
general and have been used for a long time in numerous codes. Whereas 
all of these methods based on intuitive feeling give good numerical 
results, their mathematical study remains often non-existent. 

The most interesting challenge is to solve the problem using only 
the crude rank-deficient under integrated stiffness matrix, the solution 
is obtained up to within an arbitrary spurious mode, and then to elimi- 
nate these modes from the solution in a post-processing operation. 

Unfortunately, even when the stiffness matrix is rank-sufficient, 
similar oscillations are observed when underintegration is used. In 
that case, the process of the excitations of modes similar to the hour- 
glass modes is not completely understood and these modes have never been 
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mathematically studied. 

In this report, we give precise mathematical justifications 
and answers to the questions previously mentioned. The next 
Section (Section 2.2) is devoted to the proof that the Stabilization 
method is mathematically justified. Then, in Section 2.3, we pre- 
sent a method which involves solving an underintegrated and not well- 
posed problem, then in a-posteriori eliminating the unknown degree of 
freedom. The proof of the accuracy of the method is given in Section 
2.4, and its numerical aspects and results are described in Section 2.5. 
In Section 2.6, we examine the case in which spurious oscillations 
cannot be predicted from the rank-deficiency of the stiffness matrix and 
we analyze why these modes may be excited. Finally, we apply the pre- 
vious considerations to the linear elasticity problem. 

It should be noted that the method and its results cannot be 
embedded in a classical elliptic theory: Strang's elllptlcity condition 

[49] is here violated and this non-elliptic method cannot be studied by 
the classical theory of finite element methods and numerical Integration 

[11, 12, 51, 52]. Note in these references that, both Ciarlet and 
Wahlbin crucially suppose the exactitude of the numerical scheme for 
the polynomials considered. This polynomial invariance plays a deci- 
sive role in their error estimations. We also refer to Girault [21, 

22] for his approach to the same kind of problem, non-elliptic be- 
cause of the use of partially underintegrated stiffness matrix, but 
where hourglass modes did not appear. 
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2.2 A-Priori Hourglass Control 

2.2.1 Introduction. This section is devoted to mathematical 
preliminaries to several methods consisting of adding a stabilization 
matrix to the underintegrated matrix. For clarity, we shall confine 
our attention to a simple model problem. Let Q be a regular domain 
in ]R with boundary 5Q and consider the model Neumann problem, 

(P Q ) Find u * u(x,y) such that 

-Au m f on £2 

|^-o on an 

dn 

2 

where f is an L (ft) -function satisfying 

I f dxdy - 0 (2.2) 

The questions of the existence and uniqueness of solutions to (2.1) 
(which are well-known) are taken up in Part II. 

We shall first consider a finite element approximation of (2.1) 
constructed using Q^-elements, i.e., four-node quadrilateral elements 
over which bilinear shape functions are used. Most of our notations 
and results are reproductions of those of Flanagan [18] and Belytschko 
[3, 5, 6]. Then we will attempt to extend our results to the Q^- 
elements (nine-node, biquadratic elements) and will indicate in which 
ways they differ from Belytschko’s [A]. 

The construction of finite element approximations of (2.1) involves 
the calculation of the stiffness matrix for a typical finite ele- 


I (2.1) 
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ment (2^ , which is given by the formula, 

K « VN C • VN dxdy (2.3) 

~ e ~ 

where N is a vector representing the bilinear or biquadratic shape 

functions in each element ft , 1 < e < E . 

e — — 

When (^-(respectively Q^-) elements are used to discretize the 

domain ft , K is a 4x4 matrix (resp. 9x9) and the N f s contain four 
~e - 

bilinear (resp. nine biquadratic) shape functions. We will distinguish 
exact-, full-, and under-integrations. The full integration is obtained 
using the number of Gauss integration points necessary to obtain the 
exact integration on regular square elements: 4 (resp. 9) points in our 

study. The under integration will involve the Gauss rule of lower order: 
1 (resp. 4) points. The stiffness matrix associated with a rule involv- 
ing k points will be denoted , k * 1,4,9 . 

Several authors [3, 5, 6 , 18, 30 ] proposed to add to the under inte- 
grated stiffness matrix a stabilization matrix which exhibits several 
special properties. In this section, we will prove that these proper- 
ties are indeed satisfied and that the exact stiffness matrix K can 

~e 

be computed by this method. This will be accomplished by first carry- 
ing out the integration (2.3) exactly. 

We first introduce some notations. Suppose that element is 

defined by the coordinates of its nodes (x*, y*) , 1 < I < p, p ■ 4 or 
9 . We introduce the isoparametric mapping from a master element 

n 1 . ll T 1 A ll 

^ “ |_"2’ 2j X |_~ 2 * 2j 

to ft such that 
e 
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x 


y 


l X 1 n (£,n) 
i-l 


l y 1 N (£, n ) 

i-l 


(2.4) 


where , 1 <_ I p , are the shape functions for the quadrilateral 
element on the master element. The node numbering convention is shown 
in Figure 1. 

The stiffness matrix is evaluated in (2.3) using the mapping 

(2.4) from ft to ft : 

e 


K 

-e 


I VN^ VN dxdy 


(2.5) 


where 


M 


• J s ™ T [ij |¥] ™ ■ ,Wn . 


is the Jacobian matrix of the mapping from ft to ft 


J is the inverse of its determinant, and where the gradients of the 
shape functions are derived with respect to the master element coordi- 
nates (£,n) . These matrices can be computed using (2.4) : 


N 


N 


( 2 . 6 ) 


M 


T dN T dN 

f V 


d£ 

T dN 

- <l n 


T dN 
r d n 


1 


(2.7) 
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where J Is the Jacobian of the mapping 


, dx 
det « 


Finally, we obtain the expression, 


K 

-e 


IT T T 

Axx A Ayy A 

*1 + *1 

y Ax y Ax 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


where A is the antisymmetric matrix 

A , M _ M M T (2.12' 

* dn d£ d£ dn 

T 

the Jacobian J can be expressed as y Ax • 

A study of expressed as in (2.11) and the properties of the 

matrix A will then enable us to study the effect of the underintegra- 
tion of the stiffness matrix. We will first concentrate on the 4 node- 
element and derive the exact expression of the stabilization matrix. 
Then we will discuss what form this matrix may take for the 9-node ele- 


ment . 
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2.2.2. The stabilization matrix for the bilinear element . For the 
bilinear element, the shape function vector can be written as 


where 




8 - ( 1 , - 1 , - 1 , 1 ) 
f T - (1, 1, -1, -1) 


t 1 - (1, l, l, l) 
h T - (1, -1, 1, -1) 


then the explicit form of the (4x4) matrix A is 


(2.13) 


(2.14) 


A - hs's T - ss ,T ) + |(sh T - hs T ) + ^(he’ 1 - s'h T ) 


(2.15) 


for (£,n) - (0, 0) , we obtain A^ which satisfies 


?V ■ ? T i| r . I ■ i (y 24 x 13 + i, 31 X 24 ) 

1 £>n-o 


(2.16) 


which is merely the area of the element , noted |ft e | . 

At this point, we can define precisely the matrix resulting from a 
1-point rule; this under Integra ted matrix, denoted by , is given by 


* T » T » T.T 

„(i) *oiJ £o ^ V? ?o 

‘ |!J e | + RTF 


(2.17) 


Also, if we denote B «(b^, b^) , the discrete approximations of 

the gradient VN evaluated at the integration point is given by 


h ‘ -V 


h ‘ ioi 


(2.18) 
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and 'therefore, (2 .17) takes the usual form 


»(1). 1 . T . . . T. 

*e 1?TT ( *1*1 + ~2~2^ 


The rank deficiency of can now be verified, 

Indeed, from (2.14) and (2.15) we note that 

V ■ 0 


■ 0 


and then 


K h « 0 
~e 


K t - 0 
~e ^ J 


(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 


Therefore, if we consider H and T the global hourglass and transla- 
tion, and the assembled underintegrated stiffness matrix, we have 





( 2 . 22 ) 


and this proves the rank deficiency of . Note that this "±1" 

pattern is independent of the regularity of the mesh and that h will 
take alternating values +1 and -1 at neighbor nodes as shown in Fig. 2. 

g ^ ftb 

Our goal will now be to calculate a matrix such that, if 

added to we obtain the exact stiffness matrix K given by 

~e 

(2.11). This expression does not seem easy to integrate, but the image 
of certain vectors mapped by this matrix can be easily computed using 
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orthogonality relations previously obtained (2.20) and the fact that 
~ 1 ^ “ ~ 2 ? " 

b*y « b*x « 0 (2.23) 


we obtain: 



(2.24) 


Equation (2.24) is not sufficient to compute because it gives only 

9 out of the 10 coefficients of K g (4 x 4, symmetric). It is enough 

to know X C K X , where t, x, y, and X form a set of independent 
- ~e - ~ ~ ^ ~ 

vectors. That is the case for X ■ h because 
det(x,y,t,h) * 4A ^ 0 

T 

provided the element is not singular. Then the knowledge of h K^h 
and the relations define uniquely . If we set 

h T K h - 16 F (2.25) 

- -e- 


then K is given by 

(1) — T 

K * K v ' + e YY 
-e -e 


whereas , 


again, given by (2.25), e is a scalar, and 


Y 



Q 


(2.26) 


(2.27) 
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While it is difficult to express c nicely as function of x and 
y , its exact value can be written 


- 1 
€ “ 4 


£(s T x) - n(s vT x )] 2 + C(s T y) - n(s ,T y)"l 2 
n K* + ay 43 X 12 + y i2 X 34 } + n(y 32 X 14 +y l4 X 23 ) 


d^dn 
(2 .28) 


We observe that for parallelogram elements, the denominator is constant 

and its value is the area of the domain Ifi I . In this case 

1 e 


24 m 


, 2 , 2 , 2 . 2 . 

(* 13 + x 24 + , 13 + y 24 ) 


(2.29) 


or for rectangular elements 

l 2 

J. 


i 2 + a 2 
c * — 


(2.30) 


12 £ £ 
x y 


where and are the lengths of the sides of the rectangular 

element. Also note that for such parallelogram elements y reduces to 
h . 

The expression (2.26) is often used to eliminate a-priori spurious 
modes for the kernel of K ? but the determination of £ remains a prob- 
lem. The choice £ m 0 leads to the underintegrated matrix and to the 
method to be studied in the next section. On the other hand, a cheaper 
way than the full integration of the whole matrix would be to fully 
integrate £ given by (2.28). This method would lead again to the 
full integration and is cheaper because it needs only one 4x4 integra- 
tion per element instead of 10. A more common practice is to take for 
e a simple value independent of the geometry of the element, which 



18 


is often the value obtained for a square 1/6 or sometimes any arbit- 
rary constant, as used in [5, 6). 


2.2.3. The stabilization matrix for the biquadratic element . In 

this section, we will study the effect of the underintegration of the 

9x9 stiffness matrix obtained in (2.11) with nine-node elements when 

a U- Gauss integration point rule is used. Whereas Belytschko et al. 

T 

JA] intuitively obtain another "y*Y " stabilization, we prove that this 
decomposition is not even valid for regular meshes. We then propose a 
decomposition derived on regular mesh. 

(4) 

But first we exhibit the spurious modes out of . For the bi- 

quadratic element, the shape function vector can be written as 


N - S C 


where S is the 9x9 matrix 



and 





2 2 2 
^ , n, Cn , ( i n , Cn , 



(2.31) 


(2.32) 


(2.33) 


The integration rule we are interested in involves four integration 
points (£ a , ri 0 ) » a “ 1* * • Associated with each of them, we denote 
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by and the corresponding matrix A and Jacobian. The under- 


integrated matrix is then 


„(4) i + 

K " *-> I 

~ e a*=l J a 


(2.34) 


and can also be written 


F I /u“ v“T , K a 

\ J ( h h + *2 *2 } 


a*=l a 


(2.35) 


where 


(- A a y) 

T 

(A a x) 


(2.36) 


generalizes (2.18) to a 4- point rule. 

(4) 

The rank deficiency of can now be verified. Indeed if we 

call t and h the vectors defined by 


t T - [i, i, i, i, i, i, i, i, i] 

h T - [l, 1, 1, 1, -1, -1, -1, -1, o] 


(2.37) 


we easily obtain: 


T T 

t * N “ t • S • C « 1 


T T 2 2 2 2 

h • N * h • S • £ - -4(£ + n - 12 C n ) 


and then differentiating these expressions and using (2.12) we get: 


A • t * 0 


A-h = -8 |^(1-12C 2 ) || + C(l-12n 2 ) “jJ 


the second expression vanishes when the point (£,ri) is one of the four 



20 


Gauss integration points of ft 

(5 a » n a ) “ (± — ± — ^7=- ) ; a ■ 1, 4 (2.38) 

TV 3 TV 3 

Therefore we have 


A q • t ** A^ • h “ 0 a ■ 1, 4 
and, consequently, 



(2.39) 


(2.40) 


(4) 

which proves the rank deficiency of lO Once again, we note that 

the pattern of h defined in (2.37) is independent of the geometry of 
the element and is therefore valid for a rectangular mesh as well as 
for irregular element meshes. 

The analysis of the decomposition of underintegrated and stabiliz- 
ing matrices for this element cannot be completed as completely as that 
for the 4-node element. However, Belytschko and co-workers [3] have in- 
tuitively arrived at a decomposition similar to (2.26) where y and e are 




i -p 4 b a 

4 

b* 


Y ° 

h - 

IT 

4 5 ** \ T 

1 . T _ 

- r i z z . 

Z2 

J 

(2.41) 



a«l a 

a-1 

a 




4 , r t 

T 1 



e « 

1 

100 

_ 1 . a , a , 

E , rlti • h + 

, a , a 

h • h 


(2.42) 


This decomposition does in fact satisfy several properties also satis- 
fied by the exact matrix. 


K . t - 0 


4 




21 



but for a simple square element , and its decomposition (2.26) do 

not coincide. Indeed, for this simple geometry, the calculation of 
(2.11) can be carried out explicitly and the polynomial in (£,n) ob- 
tained can be split into one part exactly integrated with A Gauss points, 
and another part of higher order that requires 9 points. This calcula- 
tion leads to the decomposition: 


r (9) = r ( 4) + n (jjr s • s* + t|t s • sj) 
xx xx e A 5 -7 -7 135 -9 -9 

k (9) - k (4) + n s • Sg + tjt s • sj) 

yy xx e A5 -8 -8 135 -9 -9 



(2.43) 


(2.44) 


and s n , s 0 , s are the 7th, 8th and 9th column vectors of S (2.32). 
-7 ~o -9 

These vectors correspond to the higher order of £ (2.33) that cannot 

be exactly integrated by a A-point rule. The form taken by the stabil- 

T 

ization matrix involves now three matrices ( s • 8^ , i = 7,8,9), is 
exact for a square element and cannot coincide with the decomposition 
found in [4]. Finally, we note that both decompositions were used in 
our a-posteriori control described in Section 2.7 on a regular mesh. 


* or also for a geometry for which the Jacobian is constant. 
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and optimal rates of convergence were only obtained with the decomposi- 
tion (2.43) . 

2.2.4. The stabilization matrix for a general heat transfer equa- 
tion. In this paragraph we the stabilization matrix for a slightly 

more complicated operator. The case of the linear elasticity operator 
is discussed later. 

Let us consider the case in which the operator is defined by 

A - £ T C 6 (2.45) 



A - 

T 

r c 

6 

where 

B - 

( JL 

''ax 

3 vT 

ay 

and 


/ c 

r \ 
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C - 



— 

\ c 

C J 



'12 

22 1 

Then the 

stiffness 

matrix associated 


K - 

T 

VN 

• C • VN dxdy 






(2.46) 


The generalization of the stabilization decomposition when elements 
are used can then be written 


(1) - T 

K - K v ' + e y • y 
-e ~e - 


(2.47) 


where 


* f?TT ? T 5 5 


(2.48) 


e * C- n e + (c . 0 + c_ - ) c + c _ ~ e 

11 XX 12 21 xy 22 yy 


(2.49) 
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and 


~xx 4 


"yy 4 


“xy 4 


n 


a 


IT T 2 

*j(£s-y - n s’ • y) d£d n 


IT T 2 

-j(£s • X - ns 1 • y) d£dn 


> (2.50) 


-j(£s T • x - r\s^ • x) (£ s^y - ns ,T • y) d£dn 


The quantities Y, B and J are the ones defined previously. Expres- 
sions similar to those given in (2.29) and (2.30) can be used to 
simplify £ . 

For a regular geometry, and corresponding to (2.29) and (2.30)> we 

have 


r _ 1 , 22 * 

xx = 24(0 ) (y 13 + y 24^ 

e 


1 , 2 . 2 , 

’yy ’ WT) <x 13 + x 2« ) 


r xy ’ 2MQT (x 13 y 13 + x 24 y 24> 


(2.51) 


As far as the 9-node element is concerned, the decomposition can 
be obtained only for regular elements. First we note that 


r (9) = r (4) 

- xy ^ xy 


(2.52) 


where the notations are similar to (2.43) and (2.44). Therefore, the 
decomposition can be written: 


K< 9 > - K< 4 ) 4. o r ( 1 T . -x 

-e - - e + C 11 (^5 s ? s ? + TTT ««> 


4 T, 

135 ?9 ?9 


+ n e C 22<B ?8 ?8 + IB -9 ?9> 


(2.53) 
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2,3 A-Posteriorl Hourglass Control 


2.3.1, Introduction and preliminaries. The basic ideas are more 
easily understood when demonstrated for the same simple model problem. 
We still focus on the model Neumann problem or its variational 

equivalent P , 

2 

Let ft be a regular (e.g. Lipschitz) domain in 1R with boundary 

2 

3ft and let f be a given L -function. Problem P^ is then, 

(Pq) Find u such that 


-Au « f in ft 

* 0 on 3ft 
9n 


(3,1) 


where the data f satisfies the compatibility condition f 

fdx - 0 (3.2) 

J ft 

Later we shall put further restrictions on ft and on f (e.g. we 
will need f G H(ft)). The kernel of the governing operator A * (-A , 
Tj^-) in (3.1) is, of course, the space of constants. Thus, whenever 
(3.2) holds, there exists a solution to (3.1) which is unique up to an 
arbitrary constant. 

To formulate a variational statement of problem P^ , we introduce 

* 

the spaces and inner products , 


2 

* The elements of V (and L (ft) / 1R) are cosets C v 3 such that u 6 JV] 
implies that u, v G H^(ft) (or L 2 (ft)) and v - u G 1R . Throughout this 
paper we frequently refer to functions v in V , meaning, of course, 
that v is a representative function in the coset [v] . 
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V - H 1 (ft) / R 

(u,v)^ * an inner product on V 
* Vu*Vvdx ; u, v € V 

h 

2 

(ffg)Q “ an inner product on L (fi)/ B 

« fgdx - —t fdx gdx 

k meas n h h 


(3.3) 


Three remarks are in order: 

i) The norm || • || 0 associated with the inner product (•,*) 0 is 

2 

the canonical norm on the quotient space L (ft) / B , 


II f II 0 “ini II f + x ll 2 

0 L (ft) 

X 6 K 

ii) According to Temam [50], there exists a constant 
ing only on f2 , such that 


(3,4) 

depend- 


llvllj >C 0 IMI 0 Vvev 


(3.5) 


iii) For all f satisfying the compatibility condition (3.2) and 
any v € V , we have 


(f ,v) Q = | fvdx 1 || f Il 0 I|v H 0 

ll v Hi 


(3.6) 


With these relations now established, we consider the variational 


statement of Pq as problem P : 
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(P) Find u G V such that 

(u,v ) 1 - (f.v) Q Vv G V (3.7) 

We can easily verify that any solution of P^ is a solution of P 
and, conversely, the solution of P satisfies the condition of P^ in 
at least, a distributional sense. Moreover, since the bilinear form 
(*, # )^ is continuous and coercive on V and since the linear form 
(f,*)^ is continuous on V if (and only if) f satisfies (3.2), the 
following result is an Immediate consequence of the Lax-Milgram Theorem: 
THEOREM I . Let f satisfy (3.2). 

Then there exists one and only one solution u G V to problem P 
and this solution depends continuously on the data f . □ 

We now consider a finite element approximation of the problem P. 

Let us now construct a finite element approximation of problem P. We 
begin by introducing a partition Q of Q into E finite elements 
so that 

E 

n - U n e 

e*l 

We shall assume that Q is such that it can be partitioned in this 
fashion into four-node quadrilateral elements over which bilinear shape 
functions are defined. Thus, if 

Ql(^e) * s P ace bilinear functions defined on 

we can introduce the finite-dimensional space 
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wherein, as usual, the label h is the mesh parameter (e.g. , 

h * max dia(ft )). The functions in are continuous and are still 

l<e<E e 

defined up to an arbitrary constant. 

Our finite-element approximation of problem P is embodied in the 
discrete problem, 


(P ) Find u h G V* 1 such that 
n 

z h h N h. 

(u , v ) 1 - (f, v ) Q 


(3.9) 


Vv h 6 v 


where, again, f satisfies condition (3.2). 

In analogy with Theorem V, we have: 

THEOREM II . Let f satisfy (3.2). Then there is one and only 

one solution to problem p in V* 1 and this solution depends 

h 


continuously on the data f 


□ 


In examining the convergence of such finite element approximations, 
we shall confine our attention throughout this section to regular mesh 
refinements. In such cases, we have the a priori asymptotic error 
estimates, 

I! u-u h |[ 1 = 0(h) , II u-u h || Q = 0(h 2 ) (3.10) 


2.3.2 The underintegrated problem . We now focus our attention 
on finite element approximations of problem P in which incomplete 
quadratures are used to evaluate the bilinear form (•j*)^ . To simpli- 
fy this study, we shall now introduce some additional assumptions: 
i) ft is the unit square. 
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n - (o.i) x (o.i) 

11) The finite elements are the squares. 


ft - (-^— i — ) x (J^i J-) 

ij ' N ’ N ; y N ’ 


1 < 1 , j < N 


ft 


ft 


l<i,j<N 


ij 


ill) The data f is L -integrable; e.g. 


f e l (Q) 


(3.11) 


In this case, we take 


h = ^ , dim V h - (N+l) 2 - 1 - 0(h" 2 ) 
N 


(3.12) 


li 2 h 

In P, we can replace f by f , its L -projection on V is 
n 


defined by 


/r n n x rc h. W h r 

(f , v ) = (f , v ). V V 6 V 


(3.13) 


For further use, we note that the projection satisfies 


IU h ll 0 i ll*ll 0 

and can be chosen such that 

[ f* 1 dx ■ 0 

J ft 


(3.14) 


(3.15) 


Now we turn to the issue of numerical integration of the stiff- 
nesses. Let I(*,*) denote a discrete inner product on C^(ft) defined 
by a numerical quadrature rule as follows: 
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V f .g) 


E 

e=l 

z w®f(C®)g(C®) 

j=l J ~J . 


(3.16) 


6 0 

Here are the quadrature weights and are the quadrature points 

for element e and G is the number of quadrature points used. 

Assuming that Gaussian quadrature is used, the choice G=4 (2x2 - Gauss 
rule) leads to an exact integration of the stiffnesses for each element: 


( h h T f h K T ir 

(u f v - I 4 (u * V ) * u K v 


(3.17) 


h h h 

for any u , v G V . Here K is the fully-integrated stiffness 
matrix and u and v are vectors of nodal degrees of freedom of u^ 
and v* 1 , respectively. 

Instead of the correct bilinear form in (3.18), we wish to consider 

an under integrated approximation to (*, # )^ in only one integra- 

* 

tion point per element is used: 

, h h. . T f h \ Tyd) 

( u , V ) 1 h c I 1 (u , V ) * u K V 


\j h h _ Tt h 
V u , v G V 


(3.18) 


Here is the underintegrated stiffness matrix. The difference be- 
tween (*,*)^ and (* t *)^ h (° n is denoted a f (*,0 and the 


corresponding stiffness matrix is K 


stab ** 


* Recall Section (2.2). 

stab — T — 

** Recall that K * cyy where e =1/6 for a rectangular mesh and 
y is given by (2.27). 
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f / ^ ^ / h h* 

a (u , v ) « (u , v - (u , v ) 1 h 


T stab 
u K v 


V h h _ 
u , v 6 V 


(3,19) 


The "under Integra ted problem", 


(P*) Find u h 6 V h such that 
n 


z h h. /r h h V/ h _ TT h 

(u * v >l,h " (f ' v >0 Vv 6 v 


(3.20) 


i8j in general , meaningless . This problem, in general, has no solution 
except for the special case in which f* 1 is orthogonal to the one- 
dimensional space of hourglass modes. 


S * (H 6 V h |(H, v h ) 1 h * 0 


V v h G V h } 


(3.21) 


A way to overcome this difficulty is to note that the under integration 
of the righthand side also leads to a rank-deficient linear form 
<•.*>, 


0,h * 


(f h , H) - 0 , V f h 6 v 1 

0,h w 

VH6S 


Note that if f satisfies (3.15) we also have 


(fh • l >0.h ■ 0 


(3.23) 


Therefore we now consider the under integrated problem P : 

h 

(P, ) Find ^ G V* such that 
h 


/ h hv f h h. 

(u , v ) lih - (f , v ) 0>h 


w h Hi 
v v 6 V 


(3.24) 
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where 


V* 


V h /H 


(3.25) 


We can now state and prove 

THEOREM III . There exists one and only one solution u* 1 to P . 

Proof : This is an immediate consequence of the Lax-Milgram theo- 

rem. Since 


(u h , u h ) 1>h = 0<£=>u h = Y x H + Y 2 
% -h 

we can consider (*,*), . as a norm on V . It is therefore coercive 

I»h 

and continuous on V** . As far as the continuity of the righthand side 

h h 

is concerned, a simple calculation shows that for any v in V we 
have 


(f 


v h ) I < 
V ; 0,h' - 


f h ll, 


h l 

v I , 


Also for any constants y^ and 


(f h , v h + Y x + Y 2 H) 0>h | - |(f", v n ) 0jh | 


h h. 


therefore 


l<fh - vh) o,hl i ll fh| lo II V ' 1 + Vj + y 2 h|| q vy. y 2 
il|f h ll 0 II v h + y 2 hIIj Vy 2 


i\H fh llo ll^lll.h 


Here we successively used (2.23), (2.22), (3,6) and the equivalence 

between the canonical norm of V* 1 and the norm ]| • || , □ 

X , n 

We have obtained a solution to the underintegrated problem . 



32 


This solution is unique in V , from a computational point of view it 
is defined up to within an arbitrary hourglass mode. We now need a 
projection to obtain a reasonable solution from any representative u 
chosen. 


2.3.3. Projection of the under integrated solution . In order to 

construct this projection, we remark that, since u* 1 is a solution of 

P and since H 6 v\ u* 1 satisfies 
h 


(u h , H) 1 - (f h , H) 0 


We wish to extend u* 1 to all of V h so that a new function u 


(3.26) 
h 


6 V n is obtained which contains an hourglass mode and which also satis- 
fies (4.8). Thus, if it is an operator from V* 1 into V* 1 , we define 


G h = ml* 1 - u 4 * + X Q H, X q 6 K (3.27) 

(Q h ,H) 1 = (f\ H) q 


This latter requirement determines ^ uniquely as 

x. — — , [< th . "> 0 - • H >il 

0 IlHlIj L 0 1J 


(3.29) 


so that u is uniquely determined as the function 

-h -h .. 

u * u + — H - 


II " II j " ll» 111 


It is instructive to consider a geometrical interpretation of our 
projection defined in (4.9). Note that the "component" of the fully 
integrated solution u* 1 orthogonal (in V h ) to H is (u h ,H), * (f^/rl)^, 
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as indicated in Fig. 3. The solutions u, of P, constitute 

n h 

the vectors generating a line "parallel to" the space H in the figure. 
The projection u ^ is then the vector defined by the orthogonal projec- 
tion of u onto this line. Indeed, by construction, 



At this point, we have established the following procedure for 
processing an underintegrated finite element approximation of problem P. 

i) Compute the underintegrated bilinear and linear forms (•,•), 

l , n 

cmi <f\-) 0 _ h 

ii) Solve problem P. for u ^ 1 

h 

Hi) Compute (u* 1 , H)^ 

iv) Construct the enhanced solution u* 1 using (3.30). 

Thus, this procedure involves the computation of an under integrated 
solution u^ to a reduced problem and its enrichment via a post- 

processing operation to obtain a new approximation u ^ . We shall now 
show that these post-processed solutions u* 1 converge to the exact 
solution u of problem P as the mesh is refined, and, remarkably, 
these approximations converge at precisely the same rate as the fully- 
integrated solution! 

Indeed we have: 

THEOREM IV : Let u > u* 1 and u^ be the solutions of P^ 

2 h 

P^ and P^ y let f be in L (ft) and satisfy (3.2). Let u be ob- 
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tained by the projection of u defined in (3.30). Then we have the 
following error estimates for s = 0 and 1 

II u h -u h || s < Cj_ h 2-s I!f || 0 (3.31) 

and 

II “““'ll* i C x h2 " S H f H 0 (3.32) 

The next section will prove this theorem. 


2.4. Convergence of the A-Posteriori Control 

2.4.1 Introduction . This section is devoted to the proof of 

Theorem VIII. The method of proof relies on the tensor properties of 

the bilinear element and of the Gauss integration rules. The problems 

and will be explicitly solved using an orthonormal basis of 

eigenvectors of (*,•)_, (*,•)-. < and (•,•)« , • Then we note that 

1 1 • h 0 , h 

2 2 

for a regular domain and mesh, f G L (Q) implies u G H ('fi'') and that 
II u - u h II x < Chllf Il 0 (4.1) 


Likewise, the Aubin-Nit sche method provides also 


II u 



< 


c'h 2 IU 


II 


0 


(4.2) 


By the triangle inequality, 

II u - J h H 1 iCh ||f || 0 + II u h - u l, || 1 C4.3) 

with a similar estimate in the || • || -norm . 

Thus, it suffices to estimate the relative error 
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h -h h 
e * u - u 


(4.4) 


2 1 

The L - and H -norms of this relative error will be explicitly calcu- 
lated and estimated. 


2.4.2. Some one-dimensional results . For reasons to be made 
clear in the next subsection, it is convenient to review briefly some 
results on one-dimensional piecewise-linear approximations on a uniform 
mesh for ■ (0,1) . Our aim here is to establish concrete relation- 
ships between various bilinear forms ( # i*)q ^ * an< * ^ #, *^1 h 

on spaces of piecewise-linear functions. 

Let D(k» a) and I denote the N+l-order matrices 


D(k,ct> 


(i. e. 


[~k a 

• 

• 

O 

o 

1 


1 

H* 

o 

• 

• 

• 

o 

o 

a 2k • • • 0 0 

• • • • • 

,I’» 

0 2 • • • 0 0 
• • • • • 

0 0 

• • 2k a 


0 0 • • • 2 0 

• 

O 

O 

1 

• • a k 

— 


1 

o 

o 

• 

• 

• 

o 

»-» 


(4.5) 

I 1 - (D( 1, 0) . Then, for a + 0 , one can show that 


det D(k, a) * (-a) N+1 det D(- £ , -1) 

- (-a) N+1 det D(- £) (4.6) 


where 


def 


D(k)=D(k, -1) (4.7) 

The values of k for which det D(k) vanishes are 
k^ ” cos , 0 £ i < N 


(4.8) 
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and the corresponding vectors (D(k^)v^ ** 0) are 

v. * {cos ^4"^} , 0 < j < N 
-l N — — 


(4,9) 


The significance of the above matrices is that in one dimension, 

12 2 
the discrete H (0,1)- , L (0,1)- and underintegrated L (0,l)-norms, 

on the space of piecewise linear C^-functions on a uniform mesh 

of N elements on (0,1) , 


vj - {v^ 6 C^(0,1) | v* 1 is linear on [eh, (e+l)h], e^O,, 


. ,N-l) 

(4.10) 


are associated with the matrices 


A Q = y D(l,Js) , A x = i D(l, -1) and A Q>h - | D(l,l) (4.11) 


respectively. In other words, 


|| v h |j 2 * v A v 

" s - -s- 


s - 0, 1 , (0,h) 


(4.12) 


where v is the vector of nodal values of v 

By using (4.6) through (4.8), one can verify that the numbers 

and 8, which render , - a.A n and A - 6. A singular are 
i ~U,n I~U ~I l~U 


“i ~ 


3(1 + cos —■) 
2(2 + cos ~) 


(4.13) 


i in 

_ 6 1 ~ cos T 

1 h 2 2 + cos ~ 
N 


(4.14) 


In particular, let <j>* * 4>*(x) , x 6 denote the piecewise 
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linear functions associated with the vectors : 


♦*( jh) * cos , 0 <_ l,j <_ N 


span { <J> ^ 0 <i<N " V 1 


(A. 15) 


Then, 


y h . i. , h ,i x 

(v > ♦ >0,h ■ “i (v • * >0 


(v h , - 8 t (v h , ♦S;, 


h _ „h 
v 6 V x 


(4.16) 

(4.17) 


Notice that the base functions <J> are orthogonal for each of the 
scalar products under consideration. 

The following remarks are in order: 

i) The denominators in (4.13) and (4.14) are non-zero, 

ii) For i * N, « 0 and the corresponding eigenfunction is the 

one-dimensional hourglass mode: 

( 1 . - 1 , 1 , - 1 , ...) 

iii) For i * 0, * 0 and the corresponding eigenfunction is 

constant. Then we have the condition (v*\l)^ * 0 as expected, 

2.4.3. Discrete norms for two dimensional meshes . The extension 
of the above results to two-dimensional rectangular meshes is straight- 
forward. Since the bilinear basis functions for V* 1 are tensor pro- 
ducts of piecewise linear functions of one variable, we can define 


^(x.y) - <{ 1 i (x)^(y) 


0 < i,j < N 


(A. 18) 
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Further, let us normalize these basis functions so that 


IU«II 0 - 1 


We can then establish the following: 

Lemma 1.1. For v* 1 G V* 1 t we have 


(v * ♦ >0.h " “i °J (V * ♦ >0 

(v h , <f. iJ ) 1 - (B ± + Bj) (v h , 4> ij ) 0 

(vh ’ * 1J) i,h * (a j 6 i + “ity (vh ‘ * lj) o 

ll ll 

Moreover, if arbitrary v 6 V is expressed in the form. 


h v .ij ^ 

v “ I v <J) 1 

0<i, j<N J 




V ij (v ' * >0 


then 


h it 2 _ 2 

v II o “ 1 v ij 

u 0<i,j<N 3 


||v h ||2 - l (B. + BJvJ 

0<i,j<N J J 


Proof : First note that 

, h ij. T / h 1 j . 

(v » ♦ 0,h I(v * (x) * (y) ) 

= a a f v h ^(x)^ (y) dxdy 

1 j J n 

- a t ^ (v h , ♦«>„ 


(4.19) 

(4.20) 

(4.21) 


(4.22) 


(4.23) 

(4.24) 
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We also have 


Finally 


f h a \ 

(v , ), 


f ^ A 1 ’/ \ J/ \ 

3 ^- 4> (x) <P (y) 


+ |p- (y) 4 » i (x) dxdy 


- 6 v h< J > 1 (x)4 )J (y) dxdy 

iJ n 

+ B v h 4> i (x)<t> J (y 

j J n 

- & ± + Bj) ( v h , 4> ij ) 0 


(y) dxdy 


, h . t r 9vh - 9v h A iJ\ 

(v • ♦ >i,h VaT* ♦ +ZT** ) 

- 6 iaj (v h ,<i> ij ) 0 + e j a 1 (v h ,d> ij ) 0 

The norms (4.23) and (4.24) are then directly obtained ^ 

In analogy with our remarks on the one-dimensional case, we observe 
that for i ■ j ■ N f <J>*^ * H , the two dimensional hourglass mode. Then 


(vh » H) i ■ ^ » H) 0 

h 


(4.25) 


• H > 0 .h " 0 


(4.26) 


Also, for i“j*0,<J>*l and the equilibrium condition (3.2) can 
be written 


0 


(4.27) 



41 


f - (f h , 4> i: *) 
ij V ’ v ; 0 


(4.28) 


2.A,A. Explicit resolution of and (P ^ 4- tt) . With the above 
results in hand, let us now return to the fully-integrated finite- 
element approximate problem P^ given in (3.9)* The solution u* 1 to 
that problem can be written 


h V 

U - Z U 4> 

0<i,j<N lJ 


“« ■ <uh - * lj) 0 


and since for in (3.9^1, 


(u\ (f. lj ) x - (B ± + &.) (u h , 4. ij ) 0 
■ <fh - ■ f ij 


(4.29) 


we have 


‘ij “ B i '+ "e f ij ; (i » j) 4 (0 * 0) 


(4.30) 


Using constructions similar to those in (A ,29) for the fully- 
integrated problem, we easily verify that the solution u* 1 to the 
underintegrated problem is representable in the form, 

^ - z u 4> iJ (4 

(i,j)*(N,N) XJ 
(iJ)^(O.O) 


a i “i 

U ii = TTp — + n ~ f ii » ^*3) ^ (0,0) and (N,N) (4.32) 

j a i e j a j-i ^ 


with 
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The cases (i,j) = (N,N) and (i,j) ■ (0,0) correspond to the arbi- 
trary hourglass mode and arbitrary constant, respectively. 

The projected approximation u* 1 defined by u* 1 * nu* 1 is con- 
structed so that projections of and u* 1 coincide; i.e. 


^ij = U ij (i,j) * ( ° ,0) 3nd (N ’ N) 
U N,N “ U N,N 


(4.33) 


2.4.5, Proof of theorem IV. Since the error function e = u 
u h is in , we use (4.29) and (4.31) to obtain 


I e * 
(l,j)l f (N,N) 13 
(i, j)l*(0,0) 


ij 


(4.34) 


where 


, h .ij v 

e « ■ <c • * >0 


.~h .ij , h .ij 

( u , ) Q ~ » 'P ) 0 


U — u 

ij ij 


Thus, from (4.30) and (4.33), 


a^cXj 


'ij KBj + a/j 


+ B. ] f ij 


(4.33) 


Then, using (4.13) and (4.14), e^. can be written as 


e u ‘ h K « f « 


(4.36) 


where 
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K 


U 


i7T 

K(cos — 


COS 



(4.37) 


and 


K(x,y) 


1 (l-bQ(l+y) 1 (2-hx)(2+y) 

4 (1+x) (l-y)+(14y) (1-x) 6 (24x)(l-y)+(2+y)(l-x) 


On the square S * J^-l,+lJ x £-1 ,+lJ /{ (-1 ,-l) , (1,1)} , K(*,*) is 
bounded and there exists a positive constant K such that 


|K(x,y) | < K V(x,y) 6 S 


(4.39) 


Therefore we have 


|k ±j I < K V(i,j) i (0,0) and (N,N) 


(4.40) 


and we can obtain using (3.13) and (4.23) 


h ii 2 


h Z 

(iJ)iK0,0) 

(i,j)^(N,N) 


2 2 
K 7 f.. 
ij 1J 


ih 4 K 2 ||f h || 2 <h‘ K 2 ||f Ho 


Also, after calculation and use of (4.24), (4.14) and (3.14), we have 


h „ 2 


= h Z 

(i,j)*(0,0) 

(i,j)*(N,N) 


<*1 + 6 j)Kl 2 f t 2 


112 h 2 K 2 ||f ||J 


2.5. Implementation and Numerical Results of the A-Posteriori Control 
For the Laplace Equation . 


In this section we first would like to indicate how the a-posteriori 
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control method is implemented, and how its time efficiency compares to 
the a-priori method. Then several numerical results will be given, 
illustrating the accuracy of the method and confirming the results 
obtained in the previous sections. 

2.5.1. Implementation of the a-posteriori method . First let us 
indicate that from a mathematical point of view the problem P* 1 is 
well-posed but computationally, the matrix obtained from this formula- 
tion is singular and the dimension of its kernel is 2. Consequently, we 
must pick two nodes, fix them a value, and solve. The first value fixes 
the constant mode, and the second one fixes the hourglass mode to be 
eliminated later. Let us fix u* 1 equal to zero at the origin and at 
the next point on the boundary (coordinates : h,0) (Figure 4. a). 
According to the error estimates (3.30) and (3.31), we may write 

U h - u* + XH + 0(h 2-s ) (5.1) 

and therefore, if we normalize H such that its nodal values are 0 or 
1, X measures precisely the value of at (h,0) (Figure 4.b), and 

approaches u(h,0) 

x - u(h,0) + 0(h°) (5.2) 

2 

But u(h,0) is 0(h ) for a smooth enough solution (u(0,0) =* 0 , 

oo 

8u/3r|(0,0) * 0) and using L -estimates [11] , 0 can be evaluated to 

2-c, e arbitrary. Finally, we have the estimate 

X = 0(h 2 G ) , c arbitrary (5.3) 


Also, the choice of H leads to 
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II H || Q - 1/3 

IIhIIj ■ M* (5.4) 

and therefore we obtain 

II XH II - 0(h 2_S ' e ) s * 0,1 (5.5) 

8 


and that proves that the post processor contribution XH can be neg- 
lected if the fixed nodes are chosen as indicated for this type of 
boundary condition. The error estimates of Theorem IV still hold up 
to within h Z . 

Unfortunately this remark has two major drawbacks: it supposes 

that u is smooth (u 6 H (Q)) and it is not valid to 9-node elements 
that will later be discussed. 

Before discussing the implementation of (3.30), we indicate that 
this projection can be simplified. Indeed, taking v* 1 * H in (4.25) 
we obtain 


(f? H) 

II -4-4 h ii 


II H||J 


< II f h ll 


II H || 


H 


— II f h H 

24 11 M 0 


(5.6) 


(f? H) 


II » || ^ 


r « llx i II 


l|H|l 0 

IlHlI, 


_h_ 


II f h ll, 


Therefore we have 


/ r h 


IlH 


h|| s ic| 


f h IL 


2-s 


0,1 


(5.8) 


and this term can be neglected without affecting the estimate of Theorem 
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VIII. The formula used in the post processor is then 

u h - ^ - XH (5.9) 

X - (u*, H) x || H II ~ 2 (5.10) 


In order to preserve the efficiency of the method provided by 
under integrat ion, one must find an efficient way to compute the para- 
meter X in (5.9). One way that suggests itself is to calculate the 
inner products of (5.10) using numerical integration. The use of a 
one point rule would be absurd and would lead to a ratio 0/0. The use 
of a 4 Gauss point rule has been numerically implemented and gives good 
results (similar to those to be presented next) but cost of this inte- 
gration is expensive, as shown in Table 3. This method is therefore 
rej ected . 

We shall now describe a more efficient method with related numer- 
ical results shown in the next subsection. This method relies on the 
fact that, for the bilinear element, the stiffness matrix can be decom- 
posed into two parts, one of which contains H in its kernel. The 
other part is such that the image of H is cheap to calculate. 

This decomposition proved in Section 2.2.2 can be written as 


, r exact ¥r under , — T 

K * K + e y *Y 

e e e ~e -e 

2 : 


(5.11) 


exact under 

where K and K respectively are the exact element stiffness 

*e =e 

matrix and its under-integrated form, e and v are obtained from 

e ie 


(2.27) and (2.28). In particular 

h*x 


He 


h - 


In. 




(5.12) 
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Even though given by (2.27) is still difficult to calculate, note 

that, during the element calculations, the vectors and b^ and the 

Jacobian \^ e \ are necessarily computed. Therefore Y e very easy 
to calculate. 

The inner product (u* 1 , H)^ may, therefore, be calculated by 

exact 

using the decomposition (2.1) of K . Introducing the nodal vec- 

tors iT and H associated with the function u* 1 and H , we have 


< m *, H). - U T K CXaX:t H - U T * It Y # Y T * H 
* 1 ~ ~ ~ ~ e e-leie 

- E E (Y* • U)(yI * H ) 

e -e -e ~e ~e 


where U and H are the values of U and H at the nodes of the 
~e ~e ~ 

element e . We note that if the values of H are +1 or -1 , the 

T 

scalar vector product y^ is always ±4 . Therefore 


H). - 4Z ± e l Y* u* 
1 e 6 i-1 e 6 


(5.13) 


(H, H). * 16Z £ 
i e e 


(5.14) 


These expressions are still exact since no approximation has been 

made on F If we suppose that the Jacobian of the element is approx- 
e 

lmately constant (true for parallelogram element), is simply ex- 

pressed as 


e " 12|?TT ( ~1 ~1 + fe 2 ~2> 


(5.15) 


The calculation of the approximate projection can be summarized in 
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the following algorithm: 

• Loop on Elements 

Calculate y , z using (2.2) and (2.6) 
~e e 

Calculate ±c (y T • U ) 
e e e 


J-Add 


{ 


\ 1 E e (Y I ' V 

X 2 + E e 


X - X^Xj 


• Loop on Nodes 

t 

*— *Node ^Node ± X 


Remark : The notations previously used are essentially those found 

in the work of Belytschko and co-workers [ 5,6] on stabilization methods. 

These methods rely on the decomposition (5.11) but the stabilization 
T 

term ey*y is a-priori added to the under-integrated matrix to prevent 
the spurious modes from the kernel of the stiffness matrix; whereas our 
control method uses the very same term a-posteriori, after solving with 
the under integrated matrix. Therefore, our method seems to be cheaper 
than the stabilization methods as summarized in Table 1. 

2.5.2. Numerical results . 

2. 5. 2. a. Regular mesh of 4-node elements . In order to illus- 
trate what has been stated, we have considered the Laplacian problem 

2 -2 

solved on a square domain partitioned into N (=h ) subdomains, for 

various values of N and we have studied the norms of the difference 
between the solution obtained with a full integration u* 1 (4 point rule 



Table 1. Operations Cost per Element for Both 

Stabilization and A-poster ior i Methods 
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and with under integration (1 point rule) u^* and u^ 1 (before and after 

i h — 

post processing). The results are shown as plot of Log|| u -u || or 

Log||u-u^|| in function of |Log h| , for s «= 0,1 . Data of vari- 
s 

ous regularities have been used: 

i) f^ is a C^-function, but not : 

r f^x.y) *= y(l-x) - yy if Y^x.y) 0 
f^(x,y) - y(|- (1-x) - y -|) if Yj (x,y) 0 

where the (^-discontinuity line is 
Yj^ (x,y)“ j (1-x) - y 

ii) f^ is a non-cont inuous function 


f 2 (x,y) - l 

if Y 1 (x,y) > 0 

f 2 (x,y) - -2 

if Y 1 (x,y) < 0 


where is the same as in i) . 

Remark : Both of these functions satisfy the compatibility 

condit ion (3.2) . 

Results obtained with the continuous function f^ are shown in 

Figure 5. When the solution has been treated by the post processor 

2 1 

(Fig. 3a.), bor both L and H norms, the representing points lie 

on lines of slope 2. This proves that whereas the estimate (1-13) is 

2 1 
optimal for the L norm (s * 0), it is not in the H norm and seems 

to be in fact better than what was expected in our study. This does 

not affect in any case the comparison with the exact solution (1.14). 

Figure 5.b shows the comparison with the crude solution u* 1 , 

obtained with two fixed nodes, and not treated by the post-processor. 
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Slopes 2 and 1 are observed and the loss (1. instead of 2. for the H^“ 

norm) corroborates the final remark of Section 3.32. 

When the function f^ is used (Figure 6), the points show os- 

2 

c illations around two lines of slope 2. (for the L -norm) and 1.65 
(for the H^-norm) proving that the estimate (3.31) still holds (Fig. 

6. a). When the solution has not been treated by the post-processor 
(Fig. 6.b), the slope 1.65 becomes 1. as expected. 

The next series of examples was intended to study the influence 
of a singularity (at the origin) for a unit square domain regularly 
partitioned. The data functions are of the form 

f a (x,y) * r a - C , a> - 2 (5.16) 

where C is a real number chosen such that the equilibrium condition 
(2.2) is satisfied. The family {f^} represents various regularities 
of data: 

f 6 H s (ft) a > s - 1 (5.17) 

a 

The result shown in Fig. 7. a is a plot of a (regularity) versus 0 

(rate of convergence of || u* 1 - u^|] ft .) . The pattern of the (a, a) 

plot seems to show a linear increase of slope 1 towards the maximum 

2 2 

value 2 reached for f G L (a* -1) for the L -norm (s e 0) . As far as 

the H^-norm (s=l) of the error is concerned, the linear increase of 

2 

slope 1 reached 1 for f G L but keeps increasing towards 2 . This 
shows that the expected error estimate 

II u h -a h ll < c h k Ilf II , s = o,i 

o m 


(5.18) 











where 


k - 1 + min (1, m) - 8 (5.19) 

is not optimal for s - 1 and m > 0 . The estimate (3.32) remains 
optimal however* 

In conclusion, these numerical results suggest that the method is 
accurate for regular meshes with a convergence rate equal to the fully 

integrated case. 

2*5.2*b* Regular mesh of 8- and 9-node elements * Since the be- 
ginning of Part II we have not discussed the under integration of the 
stiffness matrix of the 8-node elements. It is well known that this 
matrix is not rank-deficient, and the practice of the under integration 
has been widely used with good results when the mesh is regular. Since 
there is not any spurious mode, the a-posteriori control previously 
described is not needed. 

Unfortunately, the method of proof presented in Section 2*4 cannot 
be used because this element does not possess the nice tensor product 
properties on which the method relies. The only hope for a proof of 
convergence would be to obtain the result as a by-product of a result 
for the 9-node elements. 

As far as this element is concerned (9-node element), we have 
proved (Section 2.2.3) that the under integration of this element leads 
to a rank-deficient matrix; in fact, the procedure described in Section 
2.3, for the resolution of the underintegrated problem and the projec- 
tion of its solution is completely applicable to a mesh of 9-node ele- 
ments. Thus, Theorem VII is valid and the projection defined in (3.3) 
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a convergence theorem is concerned, one can establish generalizations 
of (4.16) and (4,17) to 3-node, one-dimensional elements: there exist 

0^ , Bi , 4>* and iJj* such that 


(vh » *\,h “ °i (vh * 


(v h , i^ i ) 1 - 6 i (v h , ^ i ) () 


Vv h 


Unfortunately, the basis functions $ and are different for the 

2 1 

L -under integrated and H -norms and a lemma as Lemma 1.1 cannot be ob- 
tained . 

However, in this subsection we will show numerical results 

obtained by use of the projection (3.30) for regular meshes of 9-node 

elements. Note that two types of control have been tested with similar 

T 

results: the control only involving the term in y*Y predicted by 

T 

Belytschko [4] and the complete control calculated with * i * 

8,9. (See Section 2.2.3). The results obtained with either of them are 
similar for this operator (-A) . 

For 8 and 9-node elements, the optimal rates of convergence are 
given by 


II u-uh |lg i C h k ||f || m 


s - 0,1 


(5.20) 


where 


k * 2 + rain(l,m) - s (5.21) 

3-s 1 

and the best rates of convergence 0(h ) are obtained when f 6 H (H) . 

The results obtained with functions presenting a singularity line (such 
as the functions f^ and f 2 previously defined and others) are pre- 
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sented in Table 2 (first and second lines). We obtained 1.99 and 1.74 

2 

for a discontinuous function (f G L ) f then 2.43 and 1.97 for a con- 

1 1 12 

tinuous, not C , function (f G H ), 2.95 and 1.95 for a C , not C 

2 00 

function (f G H ) and finally 4 and 3 for a C data. Therefore, the 

2 

rates 3 and 2 are reached when f is at least H or equivalent when 

4 

the solution u is in H . In this case, the convergence rate (5.1) 
does not seem to be reached. 

The second series of data involving the singularity at the origin 

(5.16) has been tested and results are shown in Fig. 6.b. The pattern 

of the (a, a) plot shows linear increases of slope 1, the predicted 

values 3 and 2 are reached for f G H^(r) according to (5.21), but the 

2 

maximum values 4 and 3 are reached for f G H (ft) . 

2.5.2.C. Irregular mesh of 4- and 9- node elements . Finally, the 
method has been tested on the quarter unit disk shown in Fig. 8 with 


f - 


a 


r 


2 

a+2 


a > - 2 


The plot (a,o) is shown in Fig. 9 and we can ooint out! 

♦ The general pattern is respected (linear increase towards a 
maximum value) 

• The maximum values 4 and 3 (9-node elements) are reduced to 
values slightly lower than 3 and 2 . 


2.6 Excitation of Spurious Modes 

The previous sections were devoted to the study of the Laplace 
equations with Neumann boundary conditions. The choice of these bound- 
ary conditions is convenient for the analysis of the hourglass instabil- 



Table 2. 


Rate of Convergence | Log ;! e ]| [v. s I Log h! (0=0,1) 

for 8- and 9-node Elements 


I REGULARITY I 

I I 

I BOUNDARY CONDITIONS I 
I AND ELEMENT I 


- ★ 

2 1 1 1 2 1 
f G L |f £ H I f G H i f6C 

t H 1 j e H 2 j e H 3 j 


I NEUMANN, 9-NODE |1.99 12.43 12.95 14.00 

I +SPECIAL PROJECTION | 1.741 1.97| 1.941 3.00 


I NEUMANN ,8-NODE EL. 11.99 12.00 12.97 14.00 

I I 1.791 1.931 1.951 3.00 


IDIRICHLET, 9-NODE EL. 12.35 12.85 12.99 13.00 

I I 1.471 1.991 1.991 2.00 

* * * * * 

IDIRICHLET, 8-NODE EL. 1 2.30 1 2. 71 |2.99 13.00 

I I 1.461 2.121 1.991 2.00 


IMIXED, 9-NODE EL. 

I 


12.00 12.00 13.63 14.00 

I 1.671 2.281 2.841 3.00 


IMIXED ,8-NODE EL. 

I 


12.00 12.00 13.84 14.00 

I 1.741 2.271 2.841 3.00 
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Figure 8. Typical Mesh on a Quarter Circle Domain. 
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ities because these modes appear explicitly in the kernel of the under- 
integrated discrete operator. When Dirichlet conditions are applied 
on a part of the boundary, even though the kernel of the under integrated 
stiffness matrix is not rank-deficient, instabilities may appear. 

In this section we would like to study the influence the boundary 
conditions have on the solution of the underintegrated problem, and 
obtain results analogous to Theorem VIII. Also we would like to explain 
how the oscillations may be excited in certain problems. The method of 
proof is similar to that presented in Section 2.4. For various boundary 
conditions, we are able to exhibit the exact eigenvalues and eigen- 
functions of the various linear and bilinear form involved. The expla- 
nation of the excitation of oscillations will result from the comparison 
of these eigenvalues. The procedure also allows us to study the under- 
integration of the operator -A+l and the control of resulting spurious 
modes. Numerical results will illustrate the theory. 

2.6.1. The under integrated problem with Dirichlet or mixed bound- 
ary conditions . This section is devoted to a generalization of the 
results obtained in Section 2.4 to the Laplacian equation with Dirichlet 
or mixed Dir ichlet-Neumann boundary conditions. Only proofs for the 
Dirichlet case will be given in this section, but their equivalent for 
the mixed case can be found in Appendix A. 

The Dirichlet case is simpler than the Neumann case because the 
hourglass mode does not belong to the new approximation space defined 
to handle the boundary condition. Therefore the stiffness matrix is no 
longer singular and can be inverted directly. In the variational for- 
mulation, similar to (3.7), the projection of the data function is not 
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necessary and the problem P is written: 


(P* 1 ) 

Find 

— h h 

u 6 V 0 

such that 




h 

(u > 

V h ) 

} l,h 

(f* 1 v* 1 ) 

U * ; 0,h 

v" 6 vj 

(6.1) 


where 

vj - {v h /v h 6 C°<n> , v h | ;!lj « (^(RIJ) , 1<1,J <», 

Remarks 

i) The fact that ^ is not singular on does not make 

the problem classically elliptic in the sense that the constant in the 
Lax Milgram Theorem is h-dependent. 

ii) When Dirichlet or mixed boundary conditions are applied. 


Ker A- * Ker A. , ■ { Ol 
1 1 , n 

Thus the post processor is not justified anymore and we may compare 
directly u^ and u* 1 . 

This comparison is carried the same way as in Section 2.4 and a 

basis of the approximation-space can be obtained. One useful 

1 2 

basis is the common eigen-basis of the matrices of the H L and 
1 2 

underintegrated H - or L - inner product. Let us consider the N-l x 
N-l matrix _ 
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The values for which det D(k) vanishes are 

k ± - cos ^ l<i<N-l (6.2) 

and the corresponding vectors (D(k^)v^, ■ 0) are 

v. * (sin } (6.3) 

1 N j"l, N-l 

Let 4. 1 - 4> i (x) , x 6 [o,l] , denote the piecewise linear function 
associated with the vector : 

<^ 1 (jh) - sin 1 £ i, j <_ N-l (6.4) 

spanU 1 }^^ - V^ 0 (6.5) 

where 

v i o * {vh 6 c ° (0 * 1) » vh(0) ■ yh<1) “ 0 

v* 1 is linear on eh, (e+l)h , 0<e<N-l} 

From this point, the remainder of the proof goes as in Section 2.4 and 
the variational problem and its underintegrated formulation can be ex- 
plicitly solved and the decomposition (4.29), (4.30) and (4.31), (4.32) 
are obtained for 1 < i,j < N-l , and we finally obtain the result for 
Dirichlet boundary conditions: 

2 

THEOREM V: Let f be a function in L (£5) . Let u be the 

solution of P : 

P : Find u 6 H* / (u,v)^ * (f,v) Q 6 H 0 (6.6) 

LqA f h be thz l 2 - pnoj&cjxon oi f onto v h and loX ^ be thz. 4o£u- 


tion P : 
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P : Find u* 1 6 v|j / (IT* 1 , h - (f h , v h ) Q h Vv h 6 (6.7) 

Then we have the ^ottouxing ewoK estimate: 

H U -^H 8 i c h2 " S II f ll 0 s - 0,1 D <6.8) 


This theorem proves that the use of the under Integra ted matrix does not 
affect the rate of convergence of the solution. The method is there- 
fore accurate and efficient. 

Data with various regularity have been tested for meshes of 4-, 8-, 
and 9-node elements, with various boundary conditions. Results are 

summarized in Tables 2 and 3* They indicate that the optimal rates of 

2 2 
convergence are obtained for f € L (Q) for the 4-node case and f 6 H 

(Q) for the 8- and 9-node case. 

2.6.2. The under integrat ion of the operator -A+l . In this sub- 
section we consider the under integrat ion of the operator associated with 
the problem 

P Q : Find u 6 H 1 (f2) such that 


f - Au + u ■ f 


l 



in ft 


3ft 


on 


(6.9) 


The usual variational formulation of P„ is 

0 


P 


Find u G 


H 1 (ft) 


such that 


(u, v) x + (u, v) 0 * (f, v) 0 , V G H X (ft) (6.10) 


The results of existence, uniqueness of solutions of P are well-known 



Table 3: Rate of Convergence jLcg 

4-node Elements (s=0,l) 


!! e H s |v.s. |Log hj for 


* * 

2 1 1 1 
f e l I f e k I f e c 

1 I T I 

t n ] t K“ ! 


REGULARITY 

OPERATOR AND 
BOUNDARY CONDITIONS 


- A 


NEUMANN 
+ PROJECTION 


1.99 12.00 12.00 

1.611 2.001 2.00 


DIRICHLET 


1.99 12.00 12.00 

1.501 1.851 2.00 


MIXED 


2.00 12.00 12.00 

1.501 1.991 1.99 


NEUMANN 

♦PROJECTION 


-4 + 1 


DIK.CELET 


MIXED 


2.00 12.00 12.00 

1.501 2.001 2.00 

— * —— — ——— 

2.00 12.00 12.00 

1.501 1.851 2.00 

* * 

2.00 12.00 12.00 

1.501 1.851 2.00 

* * 
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as are those for its discrete formulation: 


: Find u ^ G V* 1 such that 


(u\ v h ) x + (u n , v n ) 0 - (f, v n ) 0 , Vv 11 6 V r 


h h 


( 6 . 11 ) 


where V* 1 is an approximation of using bilinear elements. 

The under integration of ( • , •) ^ + ( * , •) ^ leads to the following under- 
integrated problem: 


P, : Find u* 1 G such that 
n 


/— h h h h N 

(“ , v ) 1>h + (u , v ) 0h 


where the choice of approximation space 


«• vh >o.h • Vvh e ^ 


( 6 . 12 ) 


V* 1 - V h /H 


(6.13) 


is justified by 


<v h , H) + <v\ H) 0 _ h - 0 . V v » 6 V 1 


(6.14) 


Then, the method of proof used in Section 2.3 allows us to obtain the 
existence and uniqueness of u* 1 . A projection similar to (3.30) can 
be obtained by analogy: we have 


(u\ H) x + (u\ H) 0 - (f, H) 0 (6.15) 

We therefore construct the projection as: 

u h * ttu ^ 1 = u + Xq H (6.16) 

(u h , H) x + (u\ H) 0 = (f, H) 0 


(6.17) 
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This defines uniquely X^ as 


x o- 


(f,H) 0 - (u,H) 1 - (u,H) 0 
II H || J + II H || q 


(6.18) 


Similarly to what was done in Section 2.5.1, we can use (4.25), (5.6) 
through (5.8), simplify A^ without any loss of accuracy and still use 
(5.9), (5.10) for the projection 


-h — h , 
u * u - AH 


X = tf 1 , fL) ± II H II ” 2 


(5.9) (repeated) 
(5.10) (repeated) 


The proof of the convergence of u* 1 towards u is again done by 

direct calculation of u^ and u* 1 : the explicit resolution of P^ 

and (P, + tt) leads to : 
n 


u i,j “ 1 + B i + B ij 


0 < i,j < N 


(6.19) 


and 


r - 
u 


a i°i 


i,j °i°j + a i B i + Vi 13 


0 < i,J < N 
(i,J)l‘(N,N) 


^ U NN 


NN 


( 6 . 20 ) 


These decompositions allow us to obtain u* 1 - u^ as done in Section 

2 

(2.4). Provided that f 6 L (Q) , we can obtain 


||u h -G h || s < C h z ' b '||f || 0 s * 0,1 


2-s 


( 6 . 21 ) 


Once again, the underintegration does not seem to affect the rate 
of convergence. The result can also be obtained with various boundary 
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conditions. Numerical results are summarized in Tables 5 and 6 and 
confirm the theory. 

2.6.3. Excitation of oscillations . The existence of spurious 
oscillations when under integrat ion is used is not only encountered 
when Neumann boundary conditions are applied on the whole boundary. In 
this subsection, we analyze precisely how modes that oscillate with 
wavelength of order h are excited when underintegrat ion is used, 
whereas they are damped when the integration is exact. 

For this discussion we consider the unit square 

fl - ]' 0,1 [{o,l[ (6.22) 

discretized into NxN elements. We consider the Laplace equation on ft 


- Au ■ f in ft 

u * 0 on 9ft A (x * 0} 

9u - g on 9ft/{x * 0} 

3n 


(6.23) 


For the first time we include two kinds of load: body forces and sur- 

face loads, and we will observe separately the effects of each of them. 

The eigenfunctions associated with these particular mixed boundary 
conditions are constructed as in Section 2.4. 


X 1J (x,y) - 4- i (x)tf. j (y), l<i<N (6.24) 

0<j<fo 

where ^ is defined in (4.15) (associated with Neumann boundary con- 
ditions at both ends) and is similarly defined (see Appendix A). 

These functions are defined through sine and cosine functions and there- 
fore oscillate. Among them we will distinguish '’smooth 1 * modes with 
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longer wavelengths (0(1)) from "irregular" modes with shorter wave- 
lengths (0(h)), Smooth (respectively irregular) modes correspond to 
smaller (respectively larger) values of i or j . Examples of each 
extreme are shown in Fig. 10 for N*10. 

The resolution of the fully integrated problem leads to the 


search for coefficients 


such that 


h _ ij 

U ■ E u,i X 

l<i<N J 

0<j <N 


The basis (y / is an eigenbasis for (* f *)^ an< * therefore we have 


A ij [< f . X J) 0 + ^ 8 * X ^0,3o] 


(6.26) 


where 


with 


ij 6; + Bj 


1 ,1TT 1 

6 1 ~ COs( ~N~ ~ 2 

h 2 2 + cos(-^j- - ~ 


6 1 ' cos( ffi 

h 2 2 + cos(^J-) 

N 


(6.27) 


(6.28) 


The values have been calculated exactly with these formulae 

and their values are reported in Table 4. a for N*10. The 20 highest 
values are in the shaded zone. We clearly can observe that 

i) these values range from the highest value to 1% of this value. 


ii) these values are associated with smooth modes (tensor products 
of smooth modes) • 




72 


Table 4 : 


Arrays of Ei^er. Values A , 
Table 4 a . Array A ^ 


5 . . and A 
ij 


j 

0 

i 

1 

2 

3 

4 

5 

6 7 8 

9 

10 

1 

40.45 

4.42 

1 T 54 

0;75 

0.43 

0.27 0.18 0.13 

0.10 

0.08 

1 

1 

8 i 05 

3.07 

1 .34 

0.70 

0.41 

0.26 0.17 0.12 

0.10 

0.08 

2 

1 

2 ^ 31 * 

1.58 

0.95 

O '; 57 

0.36 

0.24 0.17 0.12 

0.09 

0.08 

3 

1 

l 702 

6.85 

0.62 

0."44 

0.30 

0.21 0.15 0.11 

0.09 

0.08 

4 

I 

0755 0 . 49 - 

6 . 411 0.32 

0.24 

0.18 0.13 0.10 

0.08 

0.07 

5 

1 

0.33 

0.31 

0.27 

0.23 

0.19 

0.15 0.12 0.09 

0.08 

0.07 

6 

1 

0.21 

0.21 

0.19 

0.17 

0.14 

0.12 0.10 0.08 

0.07 

0.06 

7 

1 

0.15 

0.14 

0.14 

0.12 

0.11 

0.10 0.08 0.07 

0.06 

0.05 

8 

1 

0.11 

0.11 

0.10 

0.10 

0.09 

0.08 0.07 0.06 

0.05 

0.05 

9 

1 

0.09 

0.09 

0.08 

0.08 

0.07 

0.07 0.06 0.05 

0.05 

0.04 

10 

1 

0.08 

0.08 

0.08 

0.07 

0.07 

0.06 0.06 0.05 

0.04 

0.04 






Table 4 b . Array 



4 

i 

1 

2 

3 

4 

5 

6 7 8 

9 

10 

J 

0 

1 

40.36 

4.34 

n ;-46 

0.“67 

0*34 

0.18 0.09 0.04 

0.01 

0.00 

1 

1 

’ 7.99 

3.02 

1.27 

' 0.62 

0.33 

0.18 0.09 0.04 

0.01 

0.00 

2 

1 

2.24 

1.53 

0.90 

0.52 

0.30 

0.17 0.09 0.04 

0.01 

0.00 

3 

1 

0.94 

0.79 

0.58 

0.391 

0.25 

0.15 0.09 0.04 

0.01 

0.00 

4 

1 

0.47 

0.43 

0 . 36 | 0.28 

0.20 

0.13 0.08 0.04 

0.01 

0.00 

5 

1 

0.25 

0.24 

0.21 

0.18 

0.14 

0.11 0.07 0.04 

0.01 

0.00 

6 

1 

0.13 

0.13 

0.12 

0.11 

0.10 

0.08 0.05 0.03 

0.01 

0.00 

7 

1 

0.06 

0.06 

0.06 

0.06 

0.05 

0.05 0.04 0.03 

0.01 

0.00 

8 

1 

0.03 

0.03 

0.03 

0.03 

0.02 

0.02 0.02 0.02 

0.01 

0.00 

9 

1 

0.01 

0.01 

0.01 

0.01 

0.01 

0.01 0.01 0.01 

0.00 

0.00 

10 

1 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 0.00 0.00 

0.00 

0.00 


i 

j 

0 1 

1 

2 

Table 4 c * Array 
3 4 5 6 7 8 

9 

10 

40.45 

4.42 

1.54 0.75 

0.43 

0.27 

0.18 

0.13 

0.10 

0.08 

1 1 

8.08 

3.11 

1.36 0.71 

0.42 

0.26 

0.18 

0.13 

0.10 

0.09 

2 1 

2.32 

1.62 

0.99 0.61 

0.39 

0.26 

0.18 

0.13 

0.10 

0.09 

3 1 

1.02 

0.87 

0.67 f 6748 

0.34 

0.24 

0.18 

0.13 

0.10 

0.09 

4 1 

0.55 

0 . 5110.44 0.37 

0.29 

0.23 

0.17 

0.14 

0.11 

0.10 

5 1 

0.33 

0.32 

0.30 0.27 

0.24 

0.20 

0.17 

0.14 

0.12 

0.11 

6 1 

0.22 

0.21 

0.21 0.20 

0.19 

0.18 

0.17 

0.16 

0.14 

0.14 

7 1 

0.15 

0.15 

0.15 0.15 

0.15 

0.16 

0.17 

0.17 

0.18 

0.19 

6 1 

0.11 

0.11 

0.11 0.12 

0.13 

0.14 

0.16 

0.20 

0.26 

0.33 

9 1 

0.09 

0.09 

0.09 0.10 

0.11 

0.13 

0.16 

0.23 

0 . 42 | 0 .$ 7 | 

10 | 

0.08 

0.08 

0.09 0.09 

0.10 

0.12 

0.16 

0.251 0.57 

4 : 57 | 
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On the other hand, the eigenvalues of irregular modes are smaller and 
because of this, these modes will be damped; only smooth modes will 
contribute in (6*25). 

When the underintegration is used, and when g is zero, the solu- 
tion u is 


with 


where 


and 


~h 

u 


ij 

l£i<N J 
0<j<N 




“lj " X J, 0 


'a 


ii 


yi + a i B j 


, 3(1 ♦ cos(f - ' 

1 2(2+ cos(-y - -^j)) 

3(1 + cos(^)) 
j 2(2 + cos(^)) 


(6.29) 


(6.30) 


(6.31) 


(6.32) 


Again, the values of have been calculated exactly and they 

are reported in Table 4b* The 20 highest values are in the shaded 
zone* The comparison between Tables 7a and b shows that these 20 values 
are approximately the same and they are associated with the same smooth 
modes. In this case, irregular modes will still be damped , and one 
can predict that no oscillation will occur. 

When a load is only applied on the boundary (f ■ 0, g 4 0) , 
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is now 


where 


u ij “ A ij^ 8, X 


(6.33) 


A ij " 0j B l + a^Bj 


(6.34) 


Again the values of are reported in Table 4.c and the 20 highest 

values are in the shaded zone. Among these 20 values, three correspond 

to very irregular modes. In particular, the third value is associated 

10,10 


with x 


Therefore, we can predict a strong contribution of 

^h 


irregular inodes within the solution u , which will show oscillations. 

Finally, one could wonder if the calculation of the boundary 
integral can be calculated such that (g, X^)q jq is damped for large 
i and j . Unfortunately, no precise method has been obtained. In par- 
ticular, if the load g is a concentrated load at (x^, y^) , then 


(g. x ^o, an “ x ^ x o* V 


(6.35) 


and this value is not necessarily zero. A procedure, consisting in 
forcing (g, x*^)n * n to zero in (6.35), can be obtained in the case 
in which the load is concentrated. This procedure is discussed in a 
more general context in the next subsection. 


2.6.4. A-priori Orthogonal izat ion of the Data . Our purpose in 
this subsection is to illustrate how simple orthogonalization consider- 
ation can drastically prevent the excitation of spurious modes. In 
particular, one interprets the existence of artificial anti-hourglass 
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forces that can be applied to damp spurious oscillations* 
In order to solve the under integrated system 


A U « F 


(6,36) 


we must have the orthogonality condition 

F € (ker A) T (6.37) 

From a practical point of view, the applied load must be orthogonal to 
the spurious mode (compare with 6.35). If this is not the case, con- 
centrated spurious forces may be induced at the nodes that have been 
fixed to solve the system modwto spurious modes. 

This is clearly demonstrated in the following example where Q 
is a square partitioned into 4 or 9 4-node elements. For the Laplace 
equation, we consider two loads concentrated in two opposite points 
and satisfying the equilibrium condition as shown in Fig. 11. When 
the number of elements is even (Fig. 11. a), the system of forces are 
orthogonal to the hourglass mode H(±l-pattern indicated at the nodes) 
and, therefore, the reactions at the fixed points are zero as desired. 
Conversely, when the number of elements is odd (Fig. 11. b), the system 
of forces is not orthogonal to H and two reactions and R^, 

appear at the fixed points such that the system F^, F^, R^ and R 2 is 
orthogonal to both translation and hourglass mode. Such considerations 
can explain the peculiar rates of convergence obtained on the unit 
square discretized with NxN elements with the data 

' -6x(y 2 -l) 2 + 4 x 3 (l-3y 2 ) if x < h 

< 

■|(y 2 -l) 2 + fc(-12x 2 + 24x - 7)(l-3y 2 ) if x > *! 


f 


(6.38) 





Figure 12. a shows the rates of convergence of u* 1 towards u* 1 in the 
2 1 h 

L / IR - and H norms, where u is obtained from the under integrated 
problem and the projection, and is the solution of the fully in- 

tegrated problem. We can clearly see the good (more than optimal) 
behavior of the rates when the mesh is even; in this case, the discon- 
tinuity line (x * h) corresponds to a mesh line. However, the quality 
of the solution deteriorates when the mesh is odd or when the discon- 
tinuity line is across the mesh and coincides with the integration 
points. Note that when f^ is averaged on x * h : 

f(h, y) - «f(V*\ g) + m", y)) (6.39) 

then all the points on both plots lay on one straight line. Also 
note that the knowledge of the exact solution 


u 


3 , 2 v 2 

X (y - 1) 


if x < h 


(6. AO) 


1 -ri(-12x 2 + 24x - 7 ) ( y 2 - 1) if x > h 
lo 

The interpretation of (6.37) has been possible when Neumann 
boundary conditions are applied. When Dirichlet conditions are ap- 
plied, any right-hand side vector can produce a solution, but essen- 
tially the same type of behavior is observed (Figure 12. b). These 
results, the interpretations of (6.35) and (6.37) suggest that for any 
problem, one must pay attention to the data and make sure that it is 
orthogonal to the spurious modes. In particular, the procedure con- 
sisting in splitting a concentrated load between neighboring elements 
seems to give good results and will be discussed for elasticity prob- 


lems in the next section. 
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2.7 The Practice of Under integration In Linear Elasticity 

We devote this section to the discussion of the effects of under- 
integration in linear elasticity. We first exhibit the kernel of the 
underintegrated operator and obtain a post-processor formula similar 
to (3.30) to a-posteriori control the spurious modes. This global 
control can only be used in a very limited number of problems, but it 
suggests a local control that gives satisfactory results for all the 
examples considered. We discuss this control, its easy implementation, 
and several numerical results. 

This study is entirely qualitative -- the basis functions obtained 
in Section 2. A cannot be used to obtain eigen-functions for the elas- 
ticity operator. Both 4- and 9-node elements are discussed with an 
emphasis on the 9-node element in the numerical studies. 

2.7.1. The kernel of the discrete underintegrated linear elasti- 
city operator . We consider the linear elasticity operator defined by 

A * 8 T C B (7.1) 

where 



and C is a 3x3 symmetric matrix, 
particularize C : 

( X+2 w A 
X X+2 w 
0 0 


In the plane strain case we may 
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In order to exhibit the spurious inodes we consider the operator A 

associated with Neumann boundary conditions. In that case, the kernel 

of A consists of the usual 2-dimensional rigid body modes denoted by 

t , t and r 
-x -y 

RBM - span {t* « (J) ; t - (J) ; r - ("*)} (7. A) 

We consider the problem 

P : Find u - ( l ) 6 ["h 1 (15)1 2 /RBM such that 
~ u^ *- -* 

[ U T 6 T C 6 v dxdy - |f T • v dxdy (7.5) 

T 

where f - ( f ^ , f 2 > is a force satisfying the equilibrium conditions 
f € RBM T (7.6) 

or equivalently 


f dxdy * f. dxdy - 0 

n 1 

(f^y - f 2 x) dxdy “0 J 


(7.7) 


The existence and uniqueness of a solution for P are well known. The 
construction of finite element approximations of (7.5) involves the cal- 
culation of the (2Nx2N) stiffness matrix K g for a typical element 

ft , which is given by the formyla 
e 


K 



B T C B N dxdy 


(7.8) 


where N is a vector representing the bilinear (N*4) or biquadratic 
(N*9) shape functions in each element ft^ , 1 e <_ E . In computa- 
tional applications is evaluated using an Integration rule: 
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where, 


K 

similar to 



L 

- I 

crl 


w 


a 


(2.36), 




C B 


a 


T 


(7.9) 


(7.10) 


and w is the weight at the integration point a . Simple rank con- 
a 

siderations allow us to predict the rank of • Indeed, since 


rk(A • B) < max(rk A, rk B) N 

► 

and rk(A + B) < rk A + rk B 


we have 


rk K < 3 L 
~e — 


(7.11) 


(7.12) 


When the full integration is used, (7.12) provides no information, 

but we know that K fu11 has the correct kernel containing only rigid 
-e 

body modes. However, when under integration (L*l) is used on 4-node 
elements, we have 


rk K <3 (7.13) 

-,e — 

Therefore, the 8x8 matrix possesses at least two spurious modes. 

In fact, only two are present. Similarly, when underintegra tion 
(L=4) is used and 8- or 9-node elements, we have 


rk K <12 
~e — 


(7.14) 


This inequality predicts one spurious mode for the 16x16 matrix associ- 
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ated with 8-node elements, but when the procedure is repeated for two 
neighboring 8-node elements, the spurious modes can no longer exist in 
the global matrix [ 53 1 • We can also interpret this elimination of the 
spurious mode by noticing that neighboring element cannot share the 
mode [14]. Also note that, in spite of failing an element stability 
test [24], there are no extraneous zero energy modes in the kernel of 
the underintegrated stiffness matrix and this element is stable. 

As far as 9-node elements are concerned, the inequality (7.14) 
indicates that the 18x18 stiffness matrix has at least three spurious 
modes. In fact, there are exactly three such modes and they can be 
shared by adjacent elements. Next the modes will be explicityly de- 
scribed. 

2.7.1*a. The Spurious Modes for 4-node Elements. Let and 


H be the two hourglass vectors defined as 


H « 
-x 


( h ) 

V 



(7.15) 


where h is the hourglass nodal displacement defined in (2.14). Then 
when L*1 , we obtain from (2.18) and (2.20) 



and similarly 


B 


1 



0 


* In this section, nodal values and associated functions will be denoted 
by the same letter, the underlining differentiating them. The 

nodal values are expressed component by component. 
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Therefore, 


K (1 ^* H - K (1) * H * 0 
~e -x -e ~y 


(7.16) 


These element displacements can be put together to obtain two 


global spurious modes , also denoted by H and H and we have : 

-'■'X -y 


Ker K Un( ^ er * {span t,t y r y H,H) 
K .oc ~y f - -x -y 


(7.17) 


This defines entirely the kernel of the underintegrated matrix and the 
spurious modes for 4-node elements. 

Remark : In problems where symmetry is used for simplifications, 

the kernel of j^ un ^ er mus t respect the symmetry. If one axis of sym- 
metry (say, the x-axis) exists, then 


Ker K 


under 


span 


If the problem has two axes of symmetry (x- and y-axis) 


Ker K un< ^ er - {0) 


The spurious modes are eliminated by the symmetry conditions. 

2.7.1*b. The spurious modes for 9-node elements. Let and 


H be the two vectors defined as 

5* ’ <o> ", ’ ft 


(7.18) 


where h is the spurious mode of 9-node elements defined in (2.37). 
Using (2.36), (2.39) and (7.9), we easily get 


K^.H - K^.H = 0 
~e ~x -e ~y 


(7.19) 


Therefore H and H are two out of the three spurious modes of 
* ~x ~y 



84 


. We remark that the pattern (2.37) defining them does not depend 
on the geometry of the mesh. As far as the third spurious mode, denoted 
by 

w 

W - ( l ) (7.20) 

W 2 

is concerned, one can show that the equations defining it are 


if 


w ■ b aT • w * b „ 

Zl Z2 Z2 U 


aT 


, .aT 

"2 + h 




0 ; a-1,4 (7.21) 


or equivalently : 


yVwj - 0 
x T A°W2 - 0 


7 a-1,4 


T.a 
x A w 


T A a 


! “ Y A w 2 


(7.22) 


Note that for this system of 12 equations, we have 18 unknowns. If 
we add 5 orthogonality equations between W and t^, t^, r, and H y 
the system will define only one W (up to within a multiplicative fac- 
tor) . For a general geometry of , one cannot exhibit an explicit 

form for W ; however, when fi is a quadrilateral (strait-sided), 

e 

and when x and y are of the form 


x - (x x , x 2 , x 3 , x 4 , 55 (x x 


*5(x 3 + X 4 ) , *s(x 4 + 


we can prove that one candidate for 


»i - t y’ 


w 2 - -T x' J 


+ x 2 ) , + x 3 ) , 

x 3 ), WXj. + x 2 + x 3 + x 4 )) 
W can be written as 


(7.23) 


(7.24) 
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where 



2 0 2 

-4 -2 0 

2 4 2 

0 -2 -4 


-10 0 
110 
0 -1 -1 

0 0 1 



0 0 
0 0 
1 0 / 


(7.25) 


and x’ , y f are the vectors constructed with the first four components 
of x and £ . An example of W for a geometry satisfying (7.23) is 
shown in Figure 13 and can be constructed as follows: 

i) the displacement of a mid-side node is normal to the side, 

alternatively inwardly and outwardly oriented, with magnitude 
proportional to the length of the side, 

ii) the displacement of a corner is obtained by multiplication 
by -2 of the sum of the two displacements of the closest 
mid-side nodes. 

iii) the displacement of the centroid is zero. 

On a square, the pattern of W is well-known: 


/ - 2 , 2 , 2 , - 2 , 0 , - 1 , 0 , 1 , 0 \ 

W- (7.26) 

- \ 2 , 2 , - 2 , - 2 , - 1 , 0 , 1 . 0 . 0 / 


Contrary to 8-node elements, and because of the presence of 

and Hy , this mode can "propagate" from one element to another. For 

example, on a square mesh, if the nodal displacement vector Is W 

given by (7.26) on an element Oq , then the displacement vectors W 

+ 3H + t and W - 3H V - t on the elements to the right of 

~x ~ y ~y u 

and above allow us to construct a continuous global displacement 

also denoted W , on the mesh as shown in Figure 14. 


We finally have 




Figure 14. The Spurious Mode W : 


(a) Construction 


(b) The Mode on a 16 Element Mesh. 
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Ker R under - span (t , t , r, H , H , W} (7.27) 

-y — ** A -y — 

Remark : Similar to what we have with 4-node elements, the exis- 

tence of one axis of symmetry (say # the x-axis) reduces the kernel of 
the underintegrated stiffness matrix: 


Ker K 


under 


span 


{ Ex- 


5i* h 


+ H 3 ) 


(7.28) 


where 


!i - 3/2( !!x ■ !x } 

H 2 - -3/2(H y - t y ) 

g 3 - W + 2(t x - t y ) 


(7.29) 


J 

have been chosen such that the displacements of these modes are zero at 
the intersection of both axes for a square mesh. Contrary to the 4-node 
case, we still have a spurious mode when two axes of symmetry exist: 


Ker K unc * er ■ span (H^ + + H^} (7.30) 

This mode is shown in Figure 15. 

It is also important to point out that whereas the pattern of 

the spurious modes and H y are independent of both the geometry 

and the element, the mode W depends upon both of them. Moreover, we 

can see by construction on a square mesh that the amplitude varies 

strongly when we consider successive elements. In fact, the pattern we 

may observe is a succession of pattern H and H with increasing 

*^x -y 

amplitude. 

2.7.2. Global A-Posteriori Control in Linear Elasticity . In this 
subsection we wish to generalize (3.30) with regard to the discrete 
operators, using various kernels discussed in the previous subsection. 
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We consider the general case where 


Ker K under - RBM © spantH^ i - l,l) 


(7.31) 


where I may have the values 1, 2 or 3 , 


we obtained a control formula similar to 

_h -h • ( " ) ‘- V „ 

U " u a(H lt H x ) H 1 


We recall that for 1*1, 


(7.32) 


where the bilinear form a(»,») was obtained in the variational formu- 
lation of the initial problem. This projection satisfies: 

a(G\- H x ) - 0 (7.33) 

or, in other words, u* 1 is orthogonal to the spurious mode. We gener- 
alize this property to the elasticity problem by supposing the projection 
to be orthogonal to all the spurious modes. Therefore the control will 
consist of looking for I constants X^ (i ■ 1,1) such that 


„h — h . __ 

u * u - E X. H 

i-1,1 


a(u h , H ) - 0 for i-1, I 


(7.34) 


This leads to the system of I equations with I unknowns : 


Find \ , i ■ 1,1 such that 

E X, a(H H ) - a(u** , H ) , i-1, 1 (7.35) 

J-l.I J ~ ° ~ 

The computations involved in the control are computations of products 
of u* 1 and the spurious modes by themselves. The implementation of 
these computations is discussed in the next section. 
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2>7.2.a Implementation of the Spurious Modes Control . For the 
computation of the coefficients in (7.35), we again use the decomposition 


K full . K u„der + K > 


(7.36) 


where K un< * er satisfies 


under 

K • H - 0 


(7.37) 


3(5*, H ) - l U T • K** H 
e»l,E ~ 


(7.39) 


The expressions used for are next given for 4- or 9-node 

elements . 

2.7. 2. b. Control for 4-node E lements . For the operator defined 


'•2) 

with 



C 11 

C 12 

c n \ 


C 21 

C 22 

C 23 

(7.40) 

C 31 

C 32 

C 33 / 



we have the exact decomposition for any geometry of : 


„exact __under . 
K ■ K + 

-e 


a il VI a i2 VI 


a v*y a y*y 

*21 II 22 I I 


(7.41) 


C z +(C +C )e +C e ;C e +(C 10 +c )e 
11 xx 13 31 xy 33 yy 13 xx 12 33 xy 32 vy 


CL.e +(C 01 +C„)e +C_e ;C_e +(C„+C_) e +C 00 e 
31 xx 21 33 xy 23 yy 33 xx 32 23 xy 22 yy 


(7.42) 
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The vector y and the e f s are defined in Section 2 ((2.41) and 
(2.50)). For practival use, the expressions (2.51) are used for e . 
For linear isotropic linear material, C is given by (7.3) and 



( A+2 p) e + 

XX 


P e 


xy 


p e 


yy 


p e 


xy 


P e +(A+2P)e 
xx yy / 


(7.43) 


This expression of 
relationship: 



can be compared to the general strain-stress 



(A+2y)e + ye 

x y 


y 


xy 


y e 


xy 


ye +(A+2y)e 
x y 


(7.44) 


An algorithm similar to the one presented in Section 2.5.1 can be con- 
structed. It involves the computation of y , e and a , then the 
computation of a(H^, H^) and a(u^ . H^) , and finally the coeffi- 
cients A^ are obtained by resolution of a NxN system, N measuring 
the rank deficiency of j^ un< ^ er (n«1 or 2). 

2.7.2 .c. Control for 9-node Elements . In this subsection, devoted 
to 9-node elements, we first show why the results obtained by Belytschko 
are not sufficient to obtain a generalization of the linear elasticity 
problem, and then we propose an implementation of the control that 
leads to a stable solution converging to the exact solution with the 
optimal rate of convergence. However, for 9-node elements, we have 
not yet been able to obtain a computationally easy way to exhibit 
the third spurious mode, and the proposed results are only applica- 
ble to regular discretizations of a domain. 
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As far as the stabilization method proposed in [4] is concerned, 

algebra similar to that in Subsection 2. 7. 3. a leads to (7.41) where y 

was defined in 2.41. But, whereas the stabilization matrix constructed 

T 

with the submatrix yy eliminates H and H from the kernel of 

~x ~y 

the stiffness matrix, it does not take W into account. Indeed, we 
have 


a n M 


T T 

“21 M 


W - 0 


°12 ~‘~ T °11 ~'~ T 


Therefore this procedure cannot be used to control W . 

In order to obtain an accurate control, we have to consider a 
generalization of (2.43). Now we have 


a s *s 
11 ~i ~i 


a i2 ~i*~i 


a 21 ~i 2i 


a 22 ~i it 


• + 0 


for 7<i<9 
l<j<3 


where the vectors are defined in Section 2.2.3. Finally, using 

(2.43) and (2.52) , we have 


-e 


-e 


4Q 
e 

135 


[: 



C 31~7?7 


ll +C 33 ) f9?9 

<C 13 4<: 32 , !9!9 "I 


(C 22 +C 33 ) S9f9 ^ 

C 13f7f7 T l 

r C 33! 8 ?8 T C 32^8-8 

T J 

C 33f7f7 

L t 

C 23!s! 8 C 22!s!8 


(7.45) 


Similarly, for the 4-node case, the algorithm for the computations 



of the coefficients in (7.35) has been obtained and implemented. Numer 

T 

ical results agree with our presumptions concerning a "y # Y "-type of 
control and incline in favor of the decomposition (7.45). On a square 
domain discretized with NxN elements, we have calculated and compared 
the solutions obtained with full (exact) and underintegration for vari- 
ous boundary and symmetry conditions. Various examples considered are 

described in Figure 16. The rates of convergence were calculated by 

2 

comparing the error norms (4*0: L /RBM norm; 4*1: energy norm) ob- 
tained with N*5,6 and 7. We consistently obtained the rate 0(h^ A ) 
using a yy decomposition and 0(h ) with (7.45) for homogeneous 

materials under the action of gravity (f € C ) ; the order 3-4 being 
optimal, we may conclude that the method presented below is accurate. 

It is also efficient: for one second taken for the fully integrated 

stiffness matrix, only .61 are taken when the underintegration is used 

and only .05 seconds are taken for the control Also note that the 

1 2 

Example 4 in Figure 30 is also optimal (u € H /H ) . 

Unfortunately, this control is far from general. In particular, 
it can only be used when the exact shape of the spurious modes is known 
which is the case only when Neumann (traction) boundary conditions are 
applied on a domain discretized with regular square elements. 

2.7.3. A Local Control of the Spurious Modes 

2. 7. 3. a. Introduction . The major drawbacks of the global control 

are overcome by considering the procedure consisting of eliminating, 

element by element, the components of H , H and W and then of aver* 

x y 

aging the nodal values obtained in neighboring elements. For any ele- 
ment, we choose to do the following simplifications: 
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i) the nodal values of W are taken as if the element was 
strait-sided quadrilateral 

ii) C is diagonal* 

S t Ah 

iii) K is given by (7*41) or (7.45). 

These simplifications lead to simple calculations detailed in Appendix 
B for the 9-node element and that are easily implemented in the sub- 
routine listed in Appendix C. Note that the expression obtained for 
the control is uniquely geometric, does not depend upon the material 
properties and can be used in any linear or nonlinear problem. It only 
requires that the shape of the element not deviate too much from a 
quadrilateral . 

2.7.3.b. Numerical Results . The examples displayed in this sec- 
tion illustrate the efficiency and the accuracy of the local control 
previously described. However, we only consider linear elastic mater- 
ial in plane strain, on domains discretized with biquadratic (9-node) 
elements. Three examples are described. 

The first example is defined in Figure 17. a. We consider a square 
domain with one fixed side, under the action of pressure, gravity or a 
prescribed compressive displacement. Under any of these loads, a sing- 
ularity appears at the neighborhood of the origin. Figure 17. b shows 
how the under Integra ted solution behaves in the singularity region and 
how the control affects the results. Whereas the underintegrated solu- 
tion shows oscillations, the displacements obtained after control are 
smooth and similar to those obtained with the full integration of the 
stiffness matrix. The shear along a line AA 1 across the singularity 
region is also shown in Figure 18. Whereas undesirable oscillations 
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(a) 

Undeformed Configuration 

and Fully and Controlled- 



(b) 


Figure 17. Displacement of a Fixed Square : 

a) Definition of the Problem; b) Results. 
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are ooserved before control, the shear behaves properly after control. 

In the second example, we consider a ring under the action of an 
external pressure (Figure 19)* Here again the oscillations generated 
by the underintegration are damped when the control is applied. Only 
very slight oscillations remain, not exceeding 57*, and these can be 
easily interpreted: in the control, the expression taken for the mode 

W was obtained for quadrilateral elements. For the mesh considered, 
the elements are slightly bent and this difference explains these 
slight oscillations. The same domain (quarter ring) has also been dis- 
cretized using quadrilateral elements and the control of the underinte- 
grated solution has led to a displacement field without any oscillations 
and similar to the underintegrated displacements. Calculations of the 
stress along a radius show behavior identical to the previous example. 

The third example involves a concentrated force and illustrates 
our discussion concerning the excitation of spurious modes and the or- 
thogonaliza tion of the data. A point force is applied at a corner of 
a fixed side square discretized with a mesh refining in the neighbor- 
hood of the singularity. Strong oscillations appear in this region 
when underintegration is used, whereas the full integration solution 
is smooth (Figure 20). These oscillations show a pattern similar to 
the one used to construct the mode W (Figure 14): amplification of 

the mode 3H + t (respectively 3H r + t v ) along the x- (resp. y-) 

~x ~x ~y ~y 

direction. Therefore, according to the previous interpretations of 
(6.35) and (6.37), a way to prevent oscillations is to consider a sys- 
tem of loads similar to the load concentrated in a point A but ortho- 



100 



Figure 19, Quarter Ring Under External Pressure 

(a) Undeformed Configuration and Underintegrated 
Solution, and 

(b) Fully- and Controlled Underintegrated Solution. 




gonal to 3H + t . This is obtained by splitting the force into 3 
equal forces applied at A and its two closest nodes* Indeed, the dis 
placements obtained with this system of load only show slight oscilla 
tions. Finally, note that this control produces displacement fields 
similar to the fields obtained using full integration* 



PART III: APPLICATION TO NONLINEAR 


INCOMPRESSIBLE ELASTICITY 


3.1 Introduction 

In this concluding part, we analyze instabilities observed in dis- 
crete solutions of nonlinear problems in finite elasticity involving 
incompressible materials. We compare the behavior of the biquadratic 
(9-node) and isoparametric (8-node) elements associated with linear 
(P^) discontinuous pressures. Then we focus on 9-node elements to 
discuss the efficiency and accuracy of control of the underintegrated 
solution introduced in Part II for linear operators. 

Only Mooney-Rivlin materials are considered. They are charac- 
terized by the strain energy function 

o « C 1 (I 1 - 3) + C 2 (I 2 - 3) (3.1) 

where 1^ , i*l,3 are the principal invariants of the Cauchy-Green 
deformation tensor 

C » (1 + Vu) T ( 1 + yu) (3.2) 

where 1 is the unit tensor and Vu is the displacement gradient. 

They are: 

1 1 * tr C 

1 2 - *s(tr C 2 - (tr C) 2 ) (3.3) 

I ^ ■ det C 

The condition of incompressibility can be expressed as 
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1 3 - 1 (3.4) 

and is taken into account in a mixed formulation of the equilibrium 

problem by introducing a Lagrange multiplier P . The energy function 

Lftg 

o can be replaced by a : 

a Lag - - 3) + C 2 (I 2 - 3) + P( VTJ - 1) (3.5) 


For this choice of energy function, P is the hydro-static pressure. 
We consider the usual virtual work equilibrium equation: 


» J(J 0 ° L “ ■ iv 0 • /q 0 Po i • 6 S dv o 
• / 8 a 2 J • 6“ d > ■ « 


(3.6) 


The solution of this highly nonlinear problem is accomplished in this 
work using Newton's method. Details of the finite element method 
applied to this particular class of problems are discussed at length 
in the book of Oden [38]. A more recent account is given in Aly [1]. 


3.2 Behavior of 8- and 9-Node Elements (Full Ingegration) 

In this section, we briefly review some observations made by 
Miller [37] that are now clearly understood with the results 
obtained by Oden and Jacquotte [42, 43]. 

We consider a fixed side rectangular domain discretized with 
the refining mesh shown in Figure 21. a. We compress this domain 
Imposing a displacement on the top side. The displacement 
Increments are 5, 10, 15 and 17.5% of compression with respect to 
the original shape. 
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The solution of this problem is complicated by a stress singularity 
that occurs in the neighborhood of the corner A . Both 8- and 9-node 
displacement elements were tested associated with a linear discontin- 
uous Lagrange multiplier. The displacements show similar behaviors. 

As far as pressure is concerned, oscillations similar to the pattern 

of ker B* in [44] are observed (Figure 21. b) when 8-node elements 
h 

are used, whereas the pressure distribution is smoother with 9-node 
elements. Finally, note that the nodal averaging technique described 
in Section 1.8 has been tested by Miller [37] for a problem with simi- 
lar singularity and his conclusions corroborate ours from Part I. 

3.3 Control of the 9-Node Underintegrated Element 

In this section, we analyze how the displacements and pressures 
behave when underintegration is used. We noticed in Section 2.7 how 
the control obtained was only geometric. This is particularly useful 
for rubber-like materials and makes it very efficient when constitutive 
relations are numerically expensive to obtain. Also note that the cal- 
culation of the projected element solution (Appendix B) involves the 
knowled e of the mode W which is computed using the cu.fifie.nt deformed 
geometry of the element. We may foresee what will be one of the major 
drawbacks of the method: W is only exactly known for quadrilateral 

elements. When the element is too distorted, the approximation we do 
assuming it to be quadrilateral is too poor and the control is not 
accurate . 

Finally, we point out that the control is applied at each load 
increment. We present three examples. 
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3.3.1 Stretch of a Rubber Material Domain . The first example in- 
volves the stretch of a rubber square (Fig* 22). The domain is parti- 
tioned in 49 elements (450 degrees of freedom). The material is a 
Mooney Rivlin material with 


c x - 100 

C 2 - 20 


(3.7) 


and the incremental stretches are 25 t 50, 75 and 1007. of its undeformed 
configuration. Whereas the underintegrated solution shows slight 
oscillations of the displacement in the singularity region, the con- 
trolled solution is -smooth and similar to a fully integrated solution. 
As far as pressure is concerned, similar observations to the ones in 
Section 2.7.3.b for stresses can be made. However, the convergence is 
obtained slower: for the various displacement increment, 6, 4, 4 and 

4 (respectively 6, 5, 5 and 4) Newton iterations have been needed to 
obtain convergence with the full (resp. under) integration. Neverthe- 
less, the gain in time is almost 407* (673 sec. versus 413 sec.). 

3.3.2 Behavior in the Neighborhood of a Concentrated Force. This 
second example illustrates the ability of the orthogonaliza tion of the 
data to obtain an underintegrated solution. We consider the same 
Mooney Rivlin material (1.1, 1.7) and the problem defined in Figure 
23. a. The application of the force at only one point leads to the be- 
haviors : 


i) Slight oscillations are observed in the singularity region 
when full integration is performed (Fig. 23. b) 
ii) Uncontrollable oscillations appear when underintegration is 
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used: Figure 23. c shows the displacements after only 3 iterations 

for the first load increment. Later the solution diverges. 

As in the linear case discussed in Section 2.7.3.b, we split the 
concentrated force into three equal forces at the closest nodes; we 
observed that (Fig. 24) 

iii) When underintegration is performed, without control, the 
solution converges, but oscillations still develop (Fig. 

24. a). 

iv) The control of this solution is smooth and similar to the 
one obtained with full integration (Figs. 24. b and c). 

For this example, one supplemental iteration was needed in the third 
load increment and the gain in time also approaches 407.. 

3.3.3 Compression of a Fixed Rubber Material Domain . This final 
example reconsiders the problem used to compare the performances of 
the 8- and 9-node elements (Section 3.2). We consider two discretiza- 
tions of the domain with 25 and 49 elements, and displacement incre- 
ments of 5, 10, 15 and 17.5% with respect to the original shape. 

For the crude mesh, oscillations appear very soon when underintegration 
is performed, but the control easily corrects the solution and a dis- 
placement field close to the fully integrated field is obtained (Fig. 
25. a). But when the mesh is refined, the oscillations become more 
important and deform the element to a degree such that the control is 
not able to restore the shape of the element corner (Fig. 26. b). We 
interpret this lack of performance to the fact that the control has 
been exactly obtained for quadrilateral elements • In this case, the 


element sides curve and the element is too deformed. Also this lack 





e 25. Compress! 

(25 Eleme 
Underinte 
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of performance is observed when looking at the number of iterations 
required for the convergence of the Newton algorithm: they are 5, 5, 

5 and 5 (respectively 5, 5, 6, 7) for the fully (resp. under-) Integra 
tion in the 25 element mesh and 6, 6, and 6 (resp. 5, 6 and 13) in the 
49 element mesh. 



PART IV: CONCLUSIONS 


In this work, we study instabilities appearing in the finite 
element resolution of linear and nonlinear, compressible and 
incompressible elasticity. The study is carried both mathematically 
and numerically. Some of the principal conclusions are listed as 
follows: 

1) The use of under integrat ion in the stiffness matrix 
calculations results in rank-deficient stiffness matrix. These rank- 
deficiencies correspond to additional modes supplied to the rigid 
body modes that appropriately belong to the kernel of the operator. 

2) There is a significant class of problems in which, with 
appropriate filtering, it can be shown that an under integrated solution 
with hourglass control can yield very satisfactory answers, and 
produce a finite element method which has the same rate of convergence 
as the fully integrated method. The fact that this does indeed hold 
has been rigorously proved in this dissertation for a class of scalar 
elliptic boundary value problems, and numerically verified for a 
class of linear and nonlinear elasticity problems. 

The method developed in this work seems to give satisfactory 
results in a broad class of problems. Several questions have, 
however, arisen: 

• When an element is too distorted, the control cannot restore 
a reasonable shape. The accuracy of the control relies on the 
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approximation made on the mode W, supposing that the element remains 
quadrilateral. Does an exact computation of W (7.21) for very 
distorted elements give a better answer? 

* The lack of stiffness seems to slow the speed of convergence 
of the Newton 1 s scheme. Can the method be modified in order to 
increase the speed of convergence? 

• Finally, do the results generalize to three dimensional 
elasticity? 

The answer to these questions may provide a tremendous gain in 
computation time. 



APPENDIX A 


As far as mixed boundary conditions are concerned, we suppose 
that a Dirichlet boundary condition is applied at 0 and a Neumann 
boundary condition at 1 . For the interval Qo,l3 > we consider the 
NxN matrix 



The values for which det K(k) vanishes are: 

k t * cos (iJ + t)* 1 i 1 i n 

and the corresponding vectors (D(ki)vi * 0) are : 


{sin } 1 <_ j £ N 


The corresponding approximation space v! 1 , with basis is con- 

1 » *5 

structed as in (25]or in Section 2. A. Then, depending upon the sides 


where the various boundary conditions (D or N) are applied, tensor 

product of V* , V* or v! 1 . are to be considered. The results of 

y 1 1,0 1,*S 

Theorem II hold for the Mixed Problem. 
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APPENDIX B 


According to the simplification introduced in Section 2.7,3. a 
we have 

„ o /«y* T *v* T) + v*> T > 0 

stab m _e I T T T 

Z 135 \ 0 ; 3 (a ? ^ 7 + A 8 ^ 8 ) + 

Note that 

T T 

6 7 • h - • h * 0 

A 9 • h T - 12 


and that we can choose w * 


( w l , w 2 ) such that 


A 9* w l " V*2 


Then the control is 




with 


X i“ 4 9* U i /12 i “ 1 ’ 2 

U 7 »Wj)U 7 .u^)+U 8 »Wj)(* 8 »u^)+(4 7 .W2)(4 7 .U2)+U 8 »W2)U 8 » 

(A?*w^) 2 + (4g»Wj) 2 + U^w*) + ^g* w 2^ 2 


lift 



APPENDIX C 


SUBROUTINE PROJ(U.XY) 

C 

C THIS SUBROUTINE PROJECTS THE ELEMENT SOLUTION 

C ORTHOGONALLY W.R.T. HX,HY AND W 

C INPUT :U SOLUTION 

C XY NODAL COORDINATES (CURRENT CONFIGURATION) 

C OUTPUT :U PROJECTED SOLUTION 

C 

DIMENSION U(2,9),XY(2,9),W(2,9),S7(9),S8(9),S9(9),H(9) 
DIMENSION S7U(2) , S8U(2) ,S9U(2) ,S7W(2) ,S8W(2) 

INTEGER SIGN 

DATA S7U,S8U,S9U,S7W,S8W/10*0./ 

DATA S7/-1 . , 1., 0. ,-2. , 0., 2., 0./ 

DATA S8/-1 . ,-l . , 1., 1., 2., 0. ,-2. , 0., 0./ 

DATA S9/ 1., 1., 1., l.,-2.,-2.,-2.,-2., 4./ 

DATA H/ 1., 1., 1., 0./ 

SIGN— 1 
DO 1 K-1,2 
K1-3-K 

IF(K.EQ.2) SIGN-1 

W(K, 1 )-SIGN*(+3 ,*XY(K1 ,1 )-XY(Kl ,2)-XY(Kl ,3)-XY(Kl ,4) ) 
W(K,2)-SIGN*(+XY(K1 ,1 )-3.*XY(Kl ,2)+XY(Kl ,3)+XY(Kl ,4) ) 
W(K,3)-SIGN*(-XY(K1 ,1 )-XY(Kl ,2)+3.*XY(Kl ,3)-XY(Kl ,4) ) 
W(K,4)-SIGN*(+XY(Kl,l)+XY(Kl,2)+XY(Kl,3)-3.*XY(Kl,4)) 

W(K,5 )-SIGN*(XY(Kl ,3 )-XY(Kl ,4) ) 

W(K,6)-SIGN*(XY(K1 ,1 )-XY(Kl ,4) ) 

W(K,7)-SIGN*(XY(K1 ,1 )-XY(Kl ,2) ) 

W(K,8)-SIGN*(XY(K1 ,3)-XY(Kl ,2) ) 

1 W(K,9)-0. 

DO 2 K-1,9 
DO 2 J-1,2 

S7U(J)-S7U(J)+S7(K)*U(J,K) 

S8U(J)»S8U(J)+S8(K)*U(J,K) 

S9U(J)-S9U(J)+S9(K)*U(J,K) 

S7W(J)-S7W(J)+S7(K)*W(J,K) 

2 S8W(J)-S8W(J)+S8(K)*W( J,K) 

S9U(1)-S9U(1)/12. 

S9U(2)-S9U(2)/12. 

W1-(S7W( 1 )*S7U( 1 )+S8W( 2)*S8U( 2)+S7W( 2)*S7U( 2)+S8W( 1 )*S8U( 1 ) )/ 
. ( S7W( 1 )*S7W( 1 )+S8W( 2)*S8W( 2)+S7W( 2)*S7W( 2)+S8W( 1 )*S8W( 1 ) ) 

DO 3 K-1,9 
DO 3 J-1,2 

3 0(J,K)-U(J,K)-S9U(J)*(H(K)+l./3.)-Wl*W(J,K) 

RETURN 

END 
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